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We generalize the "renormalized" perturbation theory (RPT) formalism of Crocce and Scocci- 
marro [l|] to deal with multiple fluids in the Universe and here we present the complete calculations 
up to the one-loop level in the RPT. We apply this approach to the problem of following the non- 
linear evolution of baryon and cold dark matter (CDM) perturbations, evolving from the distinct 
sets of initial conditions, from the high redshift post-recombination Universe right through to the 
present day. In current theoretical and numerical models of structure formation, it is standard prac- 
tice to treat baryons and CDM as an effective single matter fluid - the so called dark matter only 
modeling. In this approximation, one uses a weighed sum of late time baryon and CDM transfer 
functions to set initial mass fluctuations. In this paper we explore whether this approach can be 
employed for high precision modeling of structure formation. We show that, even if we only follow 
the linear evolution, there is a large-scale scale-dependent bias between baryons and CDM for the 
currently favored WMAP5 ACDM model. This time evolving bias is significant (> 1%) until the 
present day, when it is driven towards unity through gravitational relaxation processes. Using the 
RPT formalism we test this approximation in the non-linear regime. We show that the non-linear 
CDM power spectrum in the 2-component fluid differs from that obtained from an effective mean- 
mass 1-component fluid by ~ 3% on scales of order k ~ 0.05 ZiMpc -1 at z — 10, and by ~ 0.5% at 
z — 0. However, for the case of the non-linear evolution of the baryons the situation is worse and 
we find that the power spectrum is suppressed, relative to the total matter, by ~ 15% on scales 
k ~ 0.05 h Mpc -1 at z — 10, and by ~ 3 — 5% at z = 0. Importantly, besides the suppression of the 
spectrum, the Baryonic Acoustic Oscillation (BAO) features are amplified for baryon and slightly 
damped for CDM spectra. If we compare the total matter power spectra in the 2- and 1-component 
fluid approaches, then we find excellent agreement, with deviations being < 0.5% throughout the 
evolution. Consequences: high precision modeling of the large-scale distribution of baryons in the 
Universe can not be achieved through an effective mean-mass 1-component fluid approximation; 
detection significance of BAO will be amplified in probes that study baryonic matter, relative to 
probes that study the CDM or total mass only. The CDM distribution can be modeled accurately 
at late times and the total matter at all times. This is good news for probes that are sensitive to 
the total mass, such as gravitational weak lensing as existing modeling techniques are good enough. 
Lastly, we identify an analytic approximation that greatly simplifies the evaluation of the full PT 
expressions, and it is better than < 1% over the full range of scales and times considered. 



I. INTRODUCTION 



In the current paradigm for structure formation in the Universe, it is supposed that there was an initial period of 
inflation, during which, quantum fluctuations were generated and inflated up to super-horizon scales; producing a 
near scale-invariant and near Gaussian set of primordial potential fluctuations. At the end of inflation the Universe is 
reheated, and particles and radiation are synthesized. In this hot early phase, cold thermal relics are also produced, 
these particles interact gravitationally and possibly through the weak interaction - dubbed Cold Dark Matter (CDM). 
Prior to recombination photons are coupled to electrons through Thompson scattering, and in turn electrons to 
baryonic nuclei through the Coulomb interaction. Baryon fluctuations are pressure supported on scales smaller than 
the sound horizon scale. Subsequent to recombination photons free-stream out of perturbations and baryons cool into 
the CDM potential wells. Structure formation then proceeds in a hierarchical way with small objects collapsing first 
and then merging to form larger objects. Eventually, sufficiently dense gas wells are accumulated and galaxies form. 
At late times the Universe switches from a decelerating phase of expansion to an accelerated phase. This is attributed 
to the non-zero constant energy-density of the vacuum. This paradigm has been dubbed: ACDM 0HI1- 
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FIG. 1: Baryon bias as a function of inverse spatial scale, where we have defined the baryon bias bb{k, z) = S^ n (k, z)/5^ n (k, z), 
and used Eqs (J37JI and fl38J from CTFl Solid through to dashed lines show results for redshifts z = {0, 1.0, 3.0, 5.0, 10.0, 20.0}. 

The present day energy-density budget for the ACDM model is distributed into several components, which in units 
of the critical density are: vacuum energy fl\.o ~ 0.73, matter il m .o ~ 0.27, neutrinos O^o < 10 -2 and radiation 
f2 r Q w 5 x 10~ 4 . The matter distribution can be subdivided further into contributions from CDM and baryons, with 
present day values tt c « 0.225 and Jib ~ 0.045. The detailed physics for the evolution of the radiation, CDM, baryon 
and neutrino fluctuations (8 T , S c ,5 h ,S"), from the early Universe through to recombination can be obtained by solving 
the Einstein-Boltzmann equations. These are a set of coupled non-linear partial differential equations, however while 
the fluctuations are small they may be linearized and solved [for a review see 0, 0] . At later times these fluctuations 
enter a phase of non-linear growth. Their evolution during this period can no longer be described accurately using 
the linear Einstein-Boltzmann theory and must be followed using higher order perturbation theory (hereafter PT) 
techniques Q or more directly through A^-body simulations [|[. 

The cross-over from the linear to the non-linear theory is, in general, not straightforward. Most of the theoretical 
and numerical approaches to this latter task were developed during the previous two decades, during which time the 
Standard CDM model was favored: essentially Einstein-de Sitter spacetime, f2 m = 1, with energy-density dominated 
by CDM > 97%, and baryons contributing < 3% to the total matter. For this case it is natural to assume that the 
CDM fluctuations dominate the gravitational potentials at nearly all times, with baryons playing little role in the 
formation of large-scale structure [see for example [l(J Hlj |. Thus, one simply requires a transfer function for the 
CDM distribution at some late time, to model the evolution of both components. 

In the latter part of the 1990s, cosmological tests pointed towards the ACDM paradigm, and it was recognized that 
baryons should play some role in shaping the transfer function of matter fluctuations (H, [H, EH E3 • However, rather 
than attempting to follow the CDM and baryons as separate fluids evolving from distinct sets of initial conditions, 
a simpler approximate scheme was adopted. This scheme is currently standard practice for all studies of structure 
formation in the Universe 

m he mm. it can be summarized as follows: 

1. Fix the cosmological model, specifying fi c and fib, and hence the fraction of baryons, / b . Solve for the evolution 
of all perturbed species using the linearized coupled Einstein-Boltzmann equations. Obtain transfer functions 
at z = 0, where CDM and baryons are fully relaxed: i.e. T c (fc, z = 0) » T h (k, z = 0), where T c and T b are the 
transfer functions of the CDM and baryons, respectively. 

2. Use the transfer functions to generate the linear matter power spectrum, normalized to the present day: i.e. 
Pgg(k, z = 0) » [(1 - f h )T c (k, z = 0) + f h T h (k, z = 0)] 2 Ak n w [T c (k, z = 0)} 2 Ak n , where the total matter 
fluctuation is 6 = (1 - / b )£ c + f h 5 h . 

3. Scale this power spectrum back to the initial time Zi, using the linear growth factor for the total matter 
fluctuation 5, which obtains from solving the equations of motion for a single fluid. 

4. Generate the initial CDM density field and assume that the baryons are perfect tracers of the CDM. 
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5. Evolve this effective CDM+baryon fluid under gravity using full non-linear equations of motion, either as a 
single fluid for dark matter only, or as a 2-component fluid if hydrodynamics are also to be followed. 

This effective description for the formation of structure becomes poor as the baryon fraction f h = fib/^m approaches 
~ 0(1), and for the currently favored WMAP5 cosmology / b « 0.17. In fact, in linear theory the gravitational 
relaxation between baryonic and dark matter components is only achieved at the level of < 1% by z — 0. Moreover, 
as can be seen in Fig. [1] there is a non-trivial scale-dependent bias between baryons and dark matter that exists right 
through to the present day [4!| . Failing, to take these biases into account may lead to non-negligible systematic errors 
in studies that attempt to obtain cosmological constraints from Large-Scale Structure (LSS) tests. With the next 
generation of LSS tests aiming to achieving 1% precision in measurements of the matter power spectrum |2(ij . 
it seems timely to investigate whether such modeling assumptions may impact inferences. Where we expect such 
systematic effects to play an important role is for studies based on baryonic physics such as: u sing 21cm emission 
from neutral hydrogen to trace mass density at the epoch of reionization z rc i on « 10 [HI, [H, [23], l24| : and studies of 
the Lyman alpha forest to probe the matter power spectrum [25l . l26j . 

In this, the first in a series of papers, we test whether this effective procedure can indeed be employed to follow, 
at high precision, the non-linear growth of structure formation in baryon and CDM fluids, from recombination right 
through to the present day. We do this by generalizing the matrix based "renormalized" perturbation theory (RPT) 
[50l | techniques of Crocce and Scoccimarro P, [5?], l28j to an arbitrary number of gravitating fluids. In this paper 
we content our selves with calculating the observables up to the one-loop level in the theory. In a future paper we 
consider the resummation to arbitrary orders (Somogyi & Smith 2009, in preparation). 

The late-time linear theory evolution of coupled baryon and CDM fluctuations, with sudden late time reionization 
of the inter-galactic medium was studied by [29j, |30( . Recently, [3l| explored a similar system using non-linear PT 
techniques. The work presented by these authors differs from that presented here in a number of fundamental ways. 
Firstly, these authors have assumed that, on large scales, baryons and CDM fluctuations evolve with the same initial 
conditions and that the only physical difference between the two fluids arises on small scales, where galaxy formation 
processes are able to provide some pressure support. They model this with a fixed comoving-scale Jeans criterion. 
In this work we are not so much interested in following how small-scale baryonic collapse feeds back on the mass 
distribution, but are more interested in following the very large-scale baryon and CDM fluctuations, initialized with 
their correct post-recombination spatial distributions, into the non-linear regime. Secondly, their work was developed 
within standard PT, the problems associated with this formalism are well documented 0]. Here, we develop the 
more successful RPT approach [l[. In passing, we note that a number of additional complimentary approaches and 
extensions to single-fluid RPT have recently been developed 
progress on realistic modeling of massive neutrinos and CDM 
that neutrinos play a negligible role. 

The paper is broken down as follows: In |IT]we set up the theoretical formalism: presenting the coupled equations 
of motion for the TV-component fluid. Then we show how, when written in matrix form, the equations may be solved 
at linear order. We discuss generalized initial conditions, and then discuss the specific case of baryons and CDM. In 
£|IIII we develop the non-linear PT expansion for the solutions, and show that they may be constructed to arbitrary 
order using a set of Feynman rules. In mV\ we discuss the two-point statistics of the fields, which we take as our 
lowest order observables. In $V] we give specific details of how we calculate the one-loop corrections in our theory. 
In g VII we present results for various quantities of interest. In Will we give details of an approximate RPT scheme, 
which matches the exact calculation to within < 1%, for all times of interest and on all scales where the perturbation 
theory is valid. Finally, in Willi we draw our conclusions. 




35l I36H . In addition, there has recently been 
lof . |41| . l42l |43| , for simplicity we shall assume 



II. EQUATIONS OF MOTION 



A. Standard form 



For an iV-component fluid in a uniformly expanding spacetime, the relevant equations of motion are, conservation 
of mass, momentum, and the Poisson equation. These may be written (i = 1, . . . ,N labels the i-th component) [Tj: 

' ■ ' ' +V-[(l + a i (x,T))v i (x ) T)]=Q; (1) 



<9t 
<9v. t (x, t) 
<9r 



+ W(r)vi(x, r) + (vi(x, r) • V)v;(x, r) = -V*(x, r) ; (2) 



N N 

V 2 $(x, r) = 4nGa 2 ^ p i (r)5 i (x, r) = -n m (T)H 2 (r) £ ^(x, r) . (3) 
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where r is conformal time dr = dt/a, with a the scale factor. Here we have introduced the density perturbations in 
each component as <5j(x, r) = Pi(x, T)/pi — 1 = Pi(x, T)/wip m — 1, where and pi are the density and mean density. 
Vi(x, t) is the peculiar velocity. H(t) = dlna/dr is the conformal Hubble parameter, related to the usual Hubble 
parameter, H = a/a, through H.(t) = aH. The Wi = Qi/£l m are the fractional contributions of each fluid species to 
the total matter density. They are constant in time, and normalized 



JV 



53 vh = i 



In the above we have also expressed the mean energy-density in each species in units of the critical density, 

_ pi(r) _ 8wGpi(T) 



(4) 



(5) 



and with subscript denoting present day quantities, f2^o = = To). Finally, the Hubble and density parameters 
of each species are related through the Friedmann equation, which can be written: 



W 2 (r)=a 2 TL 2 [fi m ,oa 3 + ^a,o + ^A',o« 2 ] I &K,o = 1 — &m,0 — ^A,o 
Taking the divergence of Eq. © and using Eq. © to eliminate $, we obtain 

ae<(x,r) 



(6) 



dr 



W(r)e, (x, r) + 53 [9, ( Vi )* (x, r)] (v< ), (x, r)] + (v, (x, r) • V) 9, (x, r) 



j.fe 



n m (r)7i 2 (r)53 % ^(x,r), 



(7) 



where we have set V • v(x, r) = 0(x, r) and J and k denote vector components. Next, we go to Fourier space. For a 
quantity A(x,t), we define the Fourier transform, A(k, r) and its inverse, as follows 



A(x,t) = I d 3 ke ikx i(k,r) & A(k,r) 



-^e- ikx A(x,r) 



(2tt) 3 



(8) 



Assuming Vj(x, r) to be irrotational, there exists a scalar field Uj(x, r), such that Vj(x, r) = — Vitj(x, r) and thus 
Vi(k,r) = -ikuj(k,r) , 9j(k, r) = k 2 Uj(k, t) =>■ Vj(k, r) = -ik ^-^- . 



(9) 



On inserting this relation into the equations of motion we find: 



dr 
da(k,r) 



0,(k. r) / d 3 k x d 3 k 2 <^(k - k x - k a ) ( 1 



ki-kj, 



^(ki.T^Oca.r) =0; 



(10) 



W(r)e i (k,r) + -O m (r)W 2 (r)53^^(k,r) 

i=i 

(ki+k 2 ) 2 k!.k 2 



dr ' ■-v'-^-'-' ' 2 



2k x k 2 



e i (k 1 ,r)G i (k 2 ,r) = 0. 



(11) 



d J ki d d k 2 <T(k-ki -k 2 ) 

Finally, we introduce the variable 0j(k, r) = — ©i(k, t)/H(t) and change the "time" variable to r\ — lna(r). Using 

(12) 



d dd?7 ddlna d 
= Ti—— . 



dr d?7 dr d?7 dr dr] 
and where d7i/dr = — H 2 [fl m /2 — Qa], (see Appendix |A"| . we obtain 
dSi(k,r)) 



dr] 



0<(M) = / d d kid d k 2 ^(k-ki-k a )a(k 2 ,ki)* < (ki l » 7 )fl i (ka I »7); 



(13) 
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di] 



i-^M + o A (,) 



JV 



d 3 k! d 3 k 2 ^(k - k x - k 2 )/3(ki, k 2 )^(kx, ?7)^(k 2 , 7?) , 
where in the above expressions we introduced the mode-coupling functions: 



a(ki,k 2 ) = 1 



ki-k 2 



/3(ki,k 2 ) 



(kt+k,) 2 ^-^ 
2k 1 k 2 



(14) 



(15) 



Matrix form 



Following [l[ we now rewrite the equations of motion in matrix form. To do this, we first introduce a 2JV-component 
"vector" of fields the transpose of which is given by 



C(M)= MM), OnM 



(16) 



where the index a = 1,2, ... , 2JV. With this notation, the equations of motion, Eqs (|13p and (|14[) . can be written in 
the compact form (repeated indices are summed over) 



$,tf (k,ij) + n oi tfi(k,T7) = / d 3 k! d 3 k 27 ^ ) c (k,k 1 ,k 2 )* b (k 1 ,r;)* c (k 2 ,7 / ). 



(17) 



We shall also introduce the additional compactification that repeated fc-vectors are to be integrated over, hence we 
may rewrite Eq. fl7|) as 



-» ft,*«(k, ^) + fi ab * b (k, 77) = 7 < 8 6 > c (k, k x , k 2 )* 6 (ki, r?)*e(k 2 , 77) 
The 2JV x 2JV matrix O ab appearing in Eq. IT7|) reads 



(18) 



a, 



o 



-§n m wi [i-^ + n A ] 



o 

■|o m w 2 



o 



ri m ?i; 2 [l 



7^mW 2 


















•|O m w 3 









o 







(19) 



where for the cases of JV = 1 and TV = 2 we simply take the top left 2x2 and 4x4 sub-matrices, respectively. 
The symmetrized vertex matrix may be written for general JV as, 

7S-DW-i)W)( k .k 1> k a ) = ^|M^(k- kl -k 2 ) 

ii-mm^ k ^ = %H^(k- kl -k 2 ) 

7S)( 2j )( 2j )( k 5 k l 5 k 2) = /3(ki,k 2 )^(k-k!-k 2 ), 



(20) 



where j = 1, 2, . . . , JV, and all other matrix elements are zero. For the sake of clarity, let us write explicitly the vertex 
matrix for JV = 1 and JV = 2: 



For JV = 1 the vertex matrix reads: 



T&(ki,ka) 



a(k 2) ki)/2 
a(ki,k a )/2 



7^ c (ki,k 2 ) 




/3(kx,k 2 ) 



(21) 
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• For N — 2 the vertex matrix reads: 



7 1 (8 6 ) c (ki,k 2 ) = 



7^(ki,k 2 ) 



a(k2,k x )/2 

a(ki,k 2 )/2 









a(k 2 ,ki)/2 

a(ki,k 2 )/2 



7^ c (ki,k 2 ) 





/3(ki,k 2 ) 











/3(ki,k 2 ) 



(22) 



In the above we defined 7^ to be the vertex matrix without the delta function for momentum conservation, 
i.e. 7 ^ c (k,k 1 ,k 2 ) = 7^ c (k 1 ,k 2 )^(k -ki - k 2 ). 



C. Solving for the propagator 



The main advantage of writing Eq. (|17[) in this compact form is, as was pointed out in that an implicit integral 
solution can be obtained by Laplace transforming in the variable rj. This approach, however, requires that we set the 
background cosmological model to be Einstein-de Sitter: i.e. {Cl m (rf) = 1, £Ia(v) = 0} ( we shall discuss generalizations 
to other cosmological models in ^IIDI) . This means that the matrix fl a b is independent of time, hence 



^( fl )* 6 (k, *) = 0W(k) + 7 i; ) c (k, ki,k 2 ) j ^* 6 (k, ai)* c (k, s- Sl ), 



(23) 



where (j)^ (k) = ^ a (k, f] = 0) denotes the initial conditions, set when the growth factor is one and £ (s) = sS a b + Q a b 
Multiplying by the matrix cr a b(s) and performing the inverse Laplace transform gives the formal solution 



*.(k,a) =9ab(v)$ ) 0<) + rd» ? ' 5ai ,(» ? -j 7 ')7^( k . k i k 2)*c(k 1 ,r 7 ')^(k 1 ,r7') 

Jo 

where gabiv) 1S the linear propagator given by 

gabijl) = [o- ab (s),s,rj\ , 



(24) 



(25) 



where C [f(s), s, 77] denotes the inverse Laplace transform of the function f(s) from the variable s to the variable 77. 
In Appendix IFl we prove that the linear propagator has the general form (J, k = 1, 2, . . . , N) 



3(2j-l)(2fc-l) 

ff(2i-l)(2fc) 
5(2j)(2fc) 



g (3c" - 5 + 2c~ 3 "/ 2 ) w fc + <5 i/fe , 



o"-e 3 " /2 W 



(e" - 5 + 5e"" /2 - e" 3 ^ 2 ) w k + (2 - 2e~ r '/ 2 ) ( 5 jfe , 
g (2e" - 5e~"/ 2 + 3e- 3 "/ 2 ) w k + e^ 2 S jk , 



(26) 



for 77 > 0, and g a b(v) = f° r V < due to causality, and furthermore g a b(v) ~ ¥ $ab as V ~~ * + - Clearly we can write 
the propagator as 



Sa&fa) = X/ ^ V 9ab,l , 



(27) 



where the summation runs over the set I € {1, 0, —1/2, —3/2} and the g a bj are constant matrices. Explicitly we find: 
• For N = 1, we have wi = 1 and the propagator reads: 

n/2 1 „rj r o o 1 „— 3«/2 To nl 

(28) 



1 


' 36^ 


+ 2e- 3 "/ 2 


2c'' - 2e~ 


377/2 " 


c'' 


" 3 


2 " 


e -3»)/2 


' -2 


2 " 


5 


36^ 


- 3e~ 3 "/ 2 


26^ + 3c- 


377/2 


~ "5" 


3 


2 


5 


3 


-3 
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in agreement with the usual expression where there are two solutions for the time evolution of the modes: a 
growing D + oc e v and a decaying mode D_ oc e _3,) / 2 . Thus for N = 1 onZy, ^a^o — <?af>, -1/2 = an d 



3 2 
3 2 



<?ab, -3/2 



2 -2 
-3 3 



For N = 2, the propagator takes the general form of Eq. (|2"7| , with all the (7^^ non-zero: 



(29) 



9ab,l 



9ab,-l/2 



3wi 2wi 3iV2 

3wi 2wi 3w2 

3u>i 2wi 3w2 

3u>i 2u>i 3u>2 

wi) 

-tui 

2wi 

-wi 



-2(1 



1 



2w 2 
2w 2 
2w 2 

2W2 

2w 2 

-2(1 -w 2 ) 

1 — 102 



5afc,0 



gab -3/2 



Wl 


2(1 - 


-W2 


-2w 2 














Wl 


-2w x 1 - 


w 2 2(1 


- 1V2) 














2wi 


— 2wi 2w2 


-2w 2 ~ 




3u>i 


3wi —3w2 


3w 2 




2wi 


—2wi 2w2 


— 2W2 




3wi 


3wi —3w2 


3w 2 . 





(30) 



Interestingly, we now see that there are more solutions for the time evolution of the modes, and besides the usual 
growing and decaying modes, there is an additional decaying mode g ao oc e~ v / 2 , and a static mode g a \, oc const. In 
Appendix |E] we give an explicit expression for the linear propagator for the specific case of N = 3. The interesting 
point to note in going from N = 2 to N = 3, is that, whilst the eigenvectors do change, no new eigenvalues are 
generated and this holds for general N. This explicitly follows from Eq. (|26D . 



D. Extension to general cosmological models 



An approximate treatment of the dependence of the PT and RPT kernels on cosmology was discussed by [27J, |44j , 
and we shall adopt the same approach here. The necessary changes to be made are that we use a new time variable, 

r, = log^ + (r), (31) 

where D + {t) is the linear growth factor for the appropriate cosmology Also, instead of working with the solution 
vector of the form Eq. P^|) , we work with vectors 

*r(M)= SiM/Zfo), MM), e N QL,r))/f(rjj\ . (32) 

where f(rj) = d log D + /d log a is the logarithmic growth factor for peculiar velocity fields. On performing these 
transformations Eq. (|17[) remains structurally the same, except for the elements of the matrix Q a b, which gain time 
dependence. This means that it is no longer straightforward to perform the Laplace transforms and so we do not 
obtain Eq. (|2~i|) . However, as was argued by [27|, most of the cosmological dependence is encoded in D+(rf), and so to 
a very good approximation we may simply take f2 a f,(f2 m , 57a) — * ^ab(^m = 1,S^a — 0), and keep Eq. (|24[) . We shall 
use this approximation throughout the remainder of this paper. 



E. Initial conditions 



We must now specify the initial conditions for the 2N variables. For cosmological structure formation an interesting 
case is when, for each fluid component, <5,(k, r\ = 0) and 0j(k, i] = 0) (i = 1,2,..., N) are proportional to the same 
random field [Hlj . However, the initial random field can differ from component to component. In this case we can 
write 



4 0) (k) 

where the u a , (a = 1, 2, ... , 2N) are simply numerical coefficients. 

For the case of the N = 1 component fluid, "pure growing" mode initial conditions can be obtained by setting 
u a = Ua = (1,1), and "pure decaying" mode ones by setting u a = u a 2 ^ = (2/3,-1). These are the eigenvectors of 
the matrix g a i,. In going to the more complex case of the N = 2 multi-fluid, there are, unfortunately, no such simple 



u 1 5 < ?\\ i ),u25 ( ?\k), 



U2jv-i% (k), -u 2 Ar(5^ ) (k) 



(33) 
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choices for the u a components to obtain pure "growing mode" or pure "decaying modes" at all fc-modes, unless all of 
the (k) are equivalent. To see why this is so, consider the eigenvectors of the linear propagator g a b, these are: 



,(2) 




,(3,1) 



/ W 2 



-Wi 





,(4,1) 



V 




(34) 



On dotting g a b with we find g a bU^ — e r 'u^a \ and we have a pure growing mode solution. However, to propagate 
the initial modes we require g a b dotted with 0[ O) (k). If we do this and set u a — uip then we find: 



Qb (ry)0i o) (k) = Y,g ab (vH 1] 4 0) 00 = e»v,WsW(k) = e^^(k) 



b 



(35) 



if and only if 5^ (k) = 5^(k) for all b. However, this situation is not physically interesting since, in the absence 
of any process such as pressure support, there would be no discernible difference between treating the iV-component 
fluid as an effective 1-component fluid. 

The more interesting physical case will be the situation where each fluid component evolves from a distinct set 



of initial conditions, i.e. S. 



(0) 



^ 6f, with (i^j) 



What then are the appropriate choices for the u a ? The distinct 



initial conditions that are of interest are: the dark matter, baryon, neutrino and photon, perturbations present in 
the post-recombination Universe. As described earlier the dynamics of these species are generally followed from 
the post-inflationary universe and through the recombination era using the coupled, linearized Einstein-Boltzmann 
equations [see for example |5|. In this case the fluctuations in fc-space are all proportional to the same random field, 
which is related to the primordial curvature perturbation, but with some fc-dependent transfer function that takes 
into account the physical processes that affects each species. We shall use CMBFAST [fj] for following the time evolution 
of the perturbations up to some large redshift say z — 100 as the initial condition for our non- linear computation; 
and the information on the different components is given as a set of transfer functions, which we denote Ti(fc). Thus 
for the iV-component fluid, we shall express our initial conditions, more specifically as: 



k) = [uiTi(fc), u 2 T 1 {k), u 2 N-iT N (k), u 2N T N {k)] <S (k) . 



(36) 



Importantly for the standard ACDM framework of adiabatic near scale- invariant fluctuations, all of the Ti(k) — > 1 as 

k — > 0. Hence, on large scales if we choose u a — Ua \ then </>i°^(k) is an eigenvector of the linear propagator and so 
we have pure growing mode initial conditions. On smaller scales there will, however, be some mixing of growing and 
decaying modes in the initial conditions. 



F. Application to a CDM and baryon fluid 

As a demonstration of the formalism, we now consider the specific case of the linear evolution under gravity of a two- 
component fluid of CDM and baryons. In this calculation we set the initial cosmological model to be that identified 
from the WMAP5 analysis 0: {il c = 0.2275, il h = 0.0456, A = 0.7268, h = 0.705, n s = 0.96, cr 8 = 0.812}; and 
assume negligible contribution from neutrinos. Thus, we have {w± = ft c /fl m = (1 — f h ),w 2 = J7b/^m = / b }- We use 
CMBFAST [6] to generate CDM and baryon transfer functions at the initial redshift Zj = 100: hence Ti(fc) = T c (fc, z = 
100) = T c (k) and T 2 (fc) = T b (fc, z = 100) = T b (fc). 

The implicit solution for 5 , a (k, rf), Eq. (f2~4")) . can be linearized, and the zeroth order term can be written: (k, 77) = 
5af,(^)0i O) (k) (c.f. Eq. gl])). Following our discussions in Sj] and EIIIE1 we are interested in the case where the initial 
conditions are those of the post-recombination era Universe. Prior to this time, electrons and photons are tightly 
coupled through Thompson scattering, and photons drag electrons and baryons out of the potential wells. Thus at the 
end of recombination, baryons and dark matter are spatially displaced from one another (c.f. Fig.[T]). In a somewhat 
simplified view, once the scattering processes are switched off, the evolution of the fluids can be captured by the two 
perturbations relaxing under gravity and so the evolving fields are simply some linear combination of the initial fields, 
with appropriate weight factors. 

Fig.[2]graphically demonstrates this. Here we show the ratio of the initial CDM (left panel) and baryon (right panel) 
perturbations evolving under the action of the linear propagator g a b(v)i divided through by the amplitude of the initial 
primordial density perturbation that sourced them. For large-scale growing mode initial conditions (u a — Ua ), then 
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FIG. 2: Evolution of the initial density perturbations for the CDM and baryons under the linear propagator, relative to the 
primordial density perturbation, as a function of spatial wavenumber and redshift. Top left and right panels: evolution of the 
CDM perturbations S c , and baryons, <5 b , respectively. Going from top to bottom, the blue solid through to dashed lines show 
results for redshifts z — {0,1.0,3.0,5.0,10.0,20.0}. The solid green line at the bottom shows the initial mode. Bottom left 
and right panels: ratio of the linearly evolved CDM and baryon perturbations to the initial perturbations evolved under the 
standard linear growth factor for matter perturbations. Line styles are as for above. 



we have: 

6UKv)/So(k) = M 0) (M)M>(fc) = [9ii(v)+9i2(v)]T c (k) + [ 3 i 3 (ry) + 5 i4(r/)]T b (fc) 

= (1 - / b )e" + 3/ b (l - 2c- r '/ 2 )] T c (k) + f h \c r > - 3 + 2c-"/ 2 ] T b (fc) ; (37) 

5UKv)/5o(k) = ^ 0) (k,77)/<5 (fc) = [gzi{ri)+gM\T c (k) + [g 33 (r,) + 534(77)] T h (k) 

= (1 - / b ) [c r ' - 3 + 2c-"/ 2 ! T c (fc) + \f b e r > + (1 - / b )(3 - 2 e -"/ 2 )l T h (k) .(38) 



In the limit of small and large rj these relations simplify to: 



su^v)/s (k) 

St(k,v)/So(k) 



T c (k) 

e" [(1 - / b )T c (fc) + f b T b (k)] 
T h (k) 

e „ [(l_/b) T c (fc)+/ b T b (fc) ] 



(»7<1) 
(??>1) 

(V « 1) 
(V » 1) 



(39) 
(40) 



where the rj 3> 1 relations, show that the fields evolve to an equilibrium state. 

On inspecting Fig.[2l we note that on large scales at the initial time (green solid lines) the two perturbations are 
perfectly correlated with the primordial fluctuation. On smaller scales we see that both the CDM and baryon modes 
are damped, and this is due to the Meszaros effect [see Hf| [4(|. Also the waves have an oscillating structure as a 
function of k, these are the Baryon Acoustic Oscillations (BAO) [for a discussion see [HI, [l4{. We also note that the 
BAO are almost absent in the dark matter distribution at Zi = 100, whereas for the baryons they are strongly present. 

As the system evolves linearly under gravity, we now see that on large scales (k < 0.01 /iMpc" 1 ) both solutions 
scale with time in accordance with the standard linear growth factor for the total matter perturbation on large scales: 
i.e., <5 m (k,ry) = limk^o[wid c (k, rj) + W2S h (k, -q)} — D + (rj)5 m (k,rio), which we obtain using the model of [47| with 
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FIG. 3: Graphical representation of Eg. (|45|l . The blob on the left hand side corresponds to &( n+1 ' (which clearly involves 
n + 2 initial conditions), while the two blobs on the right hand side correspond to \I/( n ~ m ) and $' m ' (involving n — m + 1 and 
m + 1 initial conditions) respectively. The trivalent vertex denotes the vertex matrix 7^. 



Q m = fi c + Jib- On smaller scales (fc > 0.01 /iMpc" ) the fluctuations grow with time, but now we note that the 
BAO features present in the baryon modes become damped, whilst those in the CDM grow. This can be seen more 
clearly by dividing each mode by D+(r])S^ (k) = D+(rj)T (k)5o(k) , for each component and this is what is plotted 

in the bottom panels of Fig.[2| Again, in the large scale limit, all Tj(fc) — > 1, and so gib{v) u b ~* D+(r]). Hence 
5f in (k,r])/[D + (ri)T c (k)5o(k)] — > 1, and likewise for the baryons. Thus we see that the amplitude of the BAO in the 
CDM are indeed increasing, and those in the baryons are decaying. 

Note that for the CDM distribution, by z — 5 the BAO are almost fully in place with < 1% changes to the 
oscillation amplitude as the distribution is evolved to z = 0. However between z — 10 and 5 there is a significant 
difference and one must be careful to account for these changes when modeling BAO at high redshift. We now draw 
attention to the important fact that the baryon distributions display more significant changes, there being an almost 
50% difference between the initial and final distributions on large-scales. As noted this will lead to significant baryon 
bias for statistical probes that are sensitive primarily to the distribution of baryons in the Universe and not the dark 
matter. In the following sections, we explore how non-linear evolution changes these results. 

III. PERTURBATION THEORY 
A. Solutions 

Following [l[ we look for power series solutions to the equations of motion, Eq. (fTT)) 

00 

* (k,»j)=X>£>(M) (41) 
j=o 

(in our notation the linear solution is as opposed to the more usual , 3>( 1 )). On inserting the above expression 
into the general solution of Eq. (|24[) . we find that 

*< 0) (M) 
*i 1} (M) 

*i 2) (M) 
*i" +1) (M) 

Fig. [3] shows a graphical representation of Eq. (|45]) . 

We may develop these solutions further by inserting into the above expressions the linear propagator g a b(fl) given 
by Eq. Q27[) . Thus, for the linear evolution of the initial waves, we have: 

*i°)(k,r,)=^e^ a M^ 0) (k). (46) 
1 



I" d V 'g ab ( V - 7/) 7 &(k, k x , k 2 )*W (k 1: r/)^ 0) (k 2 , rf) ; 
Jo 

Jo 



(42) 
(43) 

(44) 



/ d^fa-r/hSfrki.ka) ^-^(k^O^kW). 

•'° m=0 



(45) 
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Similarly, for the first order solution we have 



*W(k,»j) - /"d ?? ' 5Qb (ry-7 ? ')7^(k,k 1 ,k 2 ) 5cc ,( ?7 ')0 C '(ki).g^(^)^(k2) 5 
Jo 

= f d»/ J] e'^SaMT^Ot, ki, k 2 ) e mV '9cc>, m 4>f (k x ) £ e""W,» ^(k 2 
~' i m n 

= /V £ e^+^-^'^nEC^kr^^fc^to^^HkO^Cka) 
Jo , 



= e^m + n-*;^,^;^^ , (47) 



where in the last line we have introduced the function: 



am- jf*"- ft"- 1 '" lis 



(48) 



o 



and performed the k 2 integral using the Dirac delta function implicit in % cd (k, ki, k 2 ). 

In what follows we shall also need to make use of the 2nd order solution for the fields and after a similar procedure 
as above we find: 

^\k, V ) = 2££e i "/ 2 (m+p-?, (Z + r-p;7,) 3afcJ 7^(k 1 ,k~k 1 ) ffcc ,, m 

l,m p,q,r 

x W.pTg^k-k!-^)^^ . (49) 

In the above we have introduced a further auxiliary function, 

[l 2 e(h+h)v _ ( fl + i 2 ) e hv + h ] J [ hh ( h + ( fl ^ o, l 2 ± 0) 

" " ' [e^(Z 2 ?7 - 1) + 1] // 2 2 (*i=0,/ 2 ^0) 1 ; 

ry 2 /2 (^ = 0,^ = 0) 

B. Feynman rules for the PT solutions 

As was shown by [l[ there is a simple interpretation for the solutions ^^(k) given by Eqs (l42"|)-(|45l). The first 

equation in the hierarchy denotes an initial condition <j>a \ evolved under the linear propagator g a b. The second 
denotes the non-linear coupling between two initial waves of incoming momentum ki and k 2 through the interaction 
vertex 7„^,(k, ki, k 2 ), to produce an outgoing wave of momentum k = ki + k 2 , which is one order higher in the PT 
series and evolving under the action of g a f,. For the next equation in the series we see that the solution structure is 
now repeated, and this recurs for all higher order terms. Thus a convenient graphical representation for the solutions 
is obtained by "iterating" the graph in Fig. [3J until on the right-hand-side we are left with diagrams involving only 
trivalent vertices and no "blobs" . Associated to this graphical representation is a simple set of Feynman rules for 
obtaining the n-th order solution, \I>( n )(k, 77); these may be stated as follows: 

1. Draw all topologically distinct, connected tree diagrams containing n vertices and a directed line coming into 
the diagram from the right for each of n + 1 initial conditions, and a directed line going to the left out of the 
diagram for the final wave. Draw any number of directed internal lines running from one vertex to another, as 
required to give each vertex exactly two incoming and one outgoing attached line. It is most straightforward to 
label the direction of the lines by drawing arrows on them. 

2. Assign each vertex its own time variable, Sj, and each line its own momentum, ky, such that momentum is 
conserved at each vertex, with the momenta considered to flow in the direction of the arrows. For the single 
outgoing line, assign the momentum k and final time rj. 



3. Identify the various pieces of each diagram with the mathematical expressions as shown in Fig. [4] 
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a k b 

: 3ab(:>]-rf) 

V rf 

FIG. 4: Diagrammatic representation of the basic mathematical building blocks of the solutions. Left panel: the linear 
propagator g a b(n — T)'). Middle panel: the interaction vertex ^^.(ki, k/2). The sticks are not part of the vertex and are present 
only to show the correct number of incoming and outgoing lines. Right panel: the initial field cj$P (k) , denoted by the empty 
dot. The stick with the arrow serves only to show that an initial condition must always be connected to a diagram with a line 
going into the diagram, i.e. oriented away from the initial condition itself. 



ki + k : 



4 0, (k) 



4. Integrate over all intermediate times Sj, each between [0 , 77] . Integrate over all the independent wavevectors 
with the measure d 3 kj. (Note that the momenta of all internal lines and a single incoming line are fixed by 
momentum conservation, thus there will be n independent momenta to integrate over.) Sum over all internal 
field indices b, c, etc. 

5. Each diagram carries a factor of 2 r if there are precisely r vertices that are not symmetric with respect to 
interchanging their two incoming waves. 

6. The value of *i" } (k,r;) is given by the sum over the values of these diagrams. 

Fig. shows all diagrams up to n = 3, together with their symmetry factors calculated as in Step [5] above. As an 
example of the use of the Feynman rules, the expressions corresponding to the first two diagrams can be written, 

¥°\k, V )=g ab ( V )^ b 0) (k), (51) 
^(M) = J\s l9ab (r 1 -s 1 ) J d 3 k 1 4%(k 1 ,k-k 1 )g cc ,(s 1 )^(k 1 )g dd ,(s 1 )ct> t §\k-k l ), (52) 

which are nothing but Eqs (|46[) and (|4"T|) (after performing the Si integral in Eq. ([52")) ). Similarly, evaluating the third 
diagram, we obtain Eq. (|49[) . and so on. 




3r(0)(k l7? ) * (1) (M) *( 2 >(k,77) *( 3 )(k,j;) 

FIG. 5: Graphical representation of the topologically distinct Feynman diagrams for the perturbation theory up to third order. 
From left to right we show, <l/i n '(k), for n G {0, 1, 2, 3}, with the last two diagrams contributing to n = 3. 



IV. STATISTICS 



A. Two-point covariance matrix of density and velocity modes 



In cosmology rather than modeling ^ a (k) at a particular point, we are more interested in computing statistical 
averages of various products of the fields at different points in space. We shall represent these statistical averages 
as: (...}. For mean zero fields, we must have (^ r a (k)) = 0. Thus the lowest order statistic of interest to us is the 
two-point correlation function, or power spectrum, of ^ a (k), which can be defined: 



<* Q (k, f0¥ 6 (k', 77)) = P ab (k, v )6 D (k + k') 



(53) 
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B. Two-point covariance matrix of initial fields 



Throughout we will assume that the initial conditions are Gaussian. Consequently the statistical properties of the 
fields <?!>i 0) (k) are then completely characterized by the two-point covariance matrix: 



>i 0) (k)^ 0) (k'))^^(fe)^(k + k') 
For the case of a CDM and baryon fluid this matrix takes the explicit form: 



(54) 



u\[T c {k)] 2 Ul u 2 [T c {k)] 2 u lU3 T c (k)T h {k) Ul u 4 T c (k)T h (k) 

Ul u 2 [T c {k)] 2 u 2 2 [T c {k)] 2 u 2 u 3 T c \k)T h {k) u 2 u i T c {k)T h {k) 

Ul u 3 T c (k)T h (k) u 2 u 3 T c {k)T h {k) u 2 [T h (k)} 2 u 3 u A [T h {k)] 2 

u lUi T c (k)T b (k) u 2 u A T c {k)T h lk) u 3Ui [T h (k)] 2 ul[T h (k)} 2 



Vo(k) 



(55) 



where we defined Vo(k)5 D (k + k') = (<5o(k)(5o(k')), to be the power spectrum of primordial density fluctuations: 
i.e. Vq oc fc™, with n = 1 for the usual Harrison-Zel'dovich spectrum. The two-point correlation matrix of the 
initial fields may be expressed more compactly if we choose large-scale growing mode initial conditions, (u a — ~ 
(1, 1, 1, 1)), whereupon 



V^(k)=T a (k)T b (k)V (k) 



(56) 



with [T a (k)f = [T c (£0,T c (fc),T b (fc),T b (fc)]. 

For Gaussian initial conditions all of the multi-point correlators of the initial <fi^ (k) can be computed using Wick's 
theorem: all higher order correlators involving an odd number of fields vanish, while for an even number, there are 
(2n — 1)!! contributions corresponding to the number of different pairings of the In fields: 



^(kO...^^ 



all pair associations p pairs of 



(57) 



For later use, we now also introduce a diagrammatic representation of the initial power spectrum matrix, see Fig. [6] 

a b 
-k 



FIG. 6: Diagrammatic representation of the initial power spectrum matrix. The averaged pair of initial fields is denoted by 
the crossed dot, with the sticks and arrows serving only to show that each initial power spectrum will have two outgoing lines 
attached in any diagram. 



C. The non-linear propagator 

As the non-linear interactions are turned on, the evolution of the initial field breaks away from the simple scaling, 

^ a (k, rf) = gabiv)^ s )- T ne deviations can be thought of as, erasing the memory of the mode to its initial conditions. 
This concept can be made more concrete, through introducing the non- linear propagator [U. |27|: 

This means that we take the functional derivative of the full solution v E , a (k, rj) with respect to the initial condition 
(k'), and then perform the ensemble average. Note that the non-linear propagator G a b is defined with the 
momentum conserving delta function factored out. The expectation is taken since, again, it is the average evolution 
of the mode away from its initial conditions that we are interested in. 
In perturbation theory, the full solution is given as in Eq. I|4ip , and thus 

oo 

G ab (k, ?7 )=^5G«(k )? 7), (59) 

1=0 
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where 

M g(k,^ ( k-.o-( w ^y ). (», 

\ m ( k ) / 

Notice that in the above equation we only consider contributions to the propagator that come from terms in the 
expansion of ^I/ a (k, 77) with an odd number of initial conditions, i.e. only terms (k, if) with j an even number 
contribute. This owes to the fact that, for non-linear fields involving an even number of initial conditions, after 
performing the functional differentiation we arc left with products that involve an isolated initial field, which vanishes 
on taking the expectation: (0„ (k)) = 0. Recall that we take j — as the first term in the perturbative series for 
ty a . Thus the total non-linear propagator up to one-loop level can be written: 

G^(k,r 1 )=g ab ( V )+5G^(k,r ] ) , (61) 

where we used the fact that SG^ b '(k, rj) — g a b(v)i an d where [27l ]. 

SG^M = 4 f dr/ [ V dv"9ab'(v-v') [ d 3 k x 7^1, k - ki)gdetf ~ v") 
Jo Jo J 

x 7^(-k 1 ,k) ffcc KV)3//KV / )ff 9 6(r? / 0(ki)pW (fci) • (62) 

In Appendix [Bl we present full details of the calculation for Eq. explicitly for j = 0, 1, 2. 

As for the case of the perturbative solutions to the series expansion for the non- linear propagator with arbitrary 
numbers of loops can also be constructed in a diagrammatic fashion. Assuming Gaussian initial conditions, the 
Feynman rules for computing the n-loop correction to the propagator, SG^ (k, rj) are: 

1. Draw all tree diagrams corresponding to , E , ( 2n )(k, 77) as described in ^IIIBI Recall that these diagrams involve 
2n + 1 initial fields each. 

2. For any specific tree, drop any one initial condition, labeling the incoming line thus created with the momentum 
k (this corresponds to the functional derivative and dropping the delta function). Pair the remaining 2n initial 
fields in all possible ways into n initial power spectra, as dictated by Wick's theorem (corresponding to the 
ensemble averaging). Any specific tree then leads to (2n + l)(2n — 1)!! = (2n +1)!! n-loop diagrams, not all of 
which need to be topologically distinct. Diagrams containing a "tadpole" as a subdiagram (i.e. a subdiagram 
with a single external line) may be dropped, since they are zero, as explained below. 

3. Follow Steps HHH of ^IIIBI When assigning the lines their momenta, ensure that the two lines emanating from 
an initial power spectrum carry equal and opposite momenta. Identify the various pieces of each diagram with 
the mathematical expressions as shown in Figs0] and [5J (The number of independent momenta to integrate over 
in Step 2] will be just the number of independent loops, n.) 

4. Each diagram carries the appropriate factor of 2 r inherited form the original tree diagram. As noted above, it 
is possible that a specific loop diagram arises more than once. Of course such a diagram needs to be computed 
only once, but one must take account of its multiplicity explicitly. 

5. The value of 5G ™j} (k, v/) is given by the sum over the values of all diagrams coming from all trees. 

Fig. [7] shows all diagrams that arise in the computation of the one-loop propagator, including the symmetry factor 
computed in Step [4] above. Only the first two diagrams give non- vanishing contributions. The last diagram contains a 
tadpole (marked with the dashed box) and this gives zero contribution, since for each individual tadpole momentum 
conservation implies: a(k, — k) = /3(k, — k) = 0. Hence, diagrams containing tadpoles vanish. 

The power of the above graphical approach is that if one inspects the non-vanishing one-loop diagrams and uses 
the Feynman rules, then one can immediately write down the perturbative correction, as in Eq. (|62p . 



D. Perturbative representation of the power spectrum 

A perturbative expansion for the non-linear power spectrum may be obtained through insertion of the non-linear 
fields *i n) (c.f. Eq. 05])) into Eq. Q, whereupon 

OO 

P ab (k )77 ) = ^5P«(k,7 ? ) . (63) 
1=0 
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FIG. 7: The diagrams for the one-loop non-linear propagator G' 1 ' (k, rj). The first diagram represents Eg. (|27[) ; the second 
Eq. (|62[) ; and the last diagram contains a tadpole (marked with the dashed box) and thus evaluates to zero. 



On assuming Gaussian initial conditions, then from Wick's theorem we find that only those products of fields that 
involve an even number of initial conditions contribute to the non-linear power spectrum. Consequently we may write, 
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5P«(k, ?7 )^(k + kO=^(^ ) (k^)*f W) (k / ) r?) 

3=0 

On taking all contributions to the power spectrum matrix up to one-loop, I = 1, we have: 



Pn. 



SP 



(o) 



SP, 



(i) . 



(64) 



(65) 



lh — UJ ab r " M ab ' 

where the linear theory contribution (assuming large-scale growing mode initial conditions) is given by: 

SP^(k, V )S D (k + k') = (*W(M)¥(V,rO) 

= g aa >{r ] )T a ,{k)g w {r 1 )T b ,{k)V {k)5 D {k + k!). (66) 
We may also write the above expression using conventional matrix notation: 

P (0) (k, ?? ) = g ( v )T(k)T T (k) g T ( v )V (k). (67) 
The one-loop corrections to the power spectrum are written: 

SpW(k,r))6 D (k + k') = (^(k,^^))^^ . (68) 



Consider the sum of the first and last terms, we shall call these the "reducible" contributions and denote the sum of 
these terms: <5-P^rr ^ was shown by [27| that the non-linear propagator plays the role of a Green's function in the 



two-point sense, ($ n i 



Gac 



and furthermore that this property holds order-by-order: 



(2«) ,(o) 



SG^c (<f>c4>b }- On applying the above relations, it is straightforward to see that we must have 



<5P« R (k, n)5 D (k + k') = [(*( )(k, r,)¥ b 2 \k>, r,)) + (^(k, r,)¥ b °\k', rf)) 

= ]9ac(ri)5G® (k\ ri)Vf} (k') + SGi 1 ] (k, n)9M{r,)vf} (k)] S D (k + k') 
Specializing to growing mode initial conditions, and using momentum conservation, we find 

6P$ R (k, v ) = [gMSGQM + SG^Mgu^)] T c (k)T d (k)P (k) . 
In particular, we can write this in conventional matrix notation as 

5P^(k, n ) = [g(ry)T(fc)T T (fc)[ ( 5G( 1 )] T (k,r;)+ ( 5G( 1 )(k,^)T(fc)T T (fc)g T (r,)l V {k) . 



(69) 



(70) 



(71) 



Adding Eqs ([67f and (|7ip. we see that the sum [P*- -* (k, rf) + 5P R (k, rj)] is nothing but the one-loop truncation of 
the all-order expression 



P R (k, v ) = G(k, r?)T(fc)T T (fc)G T (k, v )V (k) . 



(72) 



The all-order reducible contribution collects all corrections to P(k,?7) that are proportional to the initial power 
spectrum at the same scale k. 
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FIG. 8: The non-zero diagrams for the one-loop power spectrum P"'(k, 77). The dashed lines indicate the points at which 
the two trees were glued together. The first diagram represents Eq. (|66[) ; the second and third graphs are the two reducible 
contributions in Eq. (|71|l ; and the fourth graph represents the mode-coupling term of Eq. (|73|l . 

Turning now to the second term in Eq. (|68p . conventionally called the "mode-coupling" contribution, we find that 
it can be written [l|: 

8P$ucQi, V) = 2 f oV f dr?" [ d 3 k x g ac {v ~ ^^(ki.k - k 1 )^(r ? ')5ee'(r/') 

Jo Jo 7 (73) 

x ff6/(r? - ^'^iC-ki.ka - k) ffsg ,(V')^'(V')^S(fe)^(|k - k x |) . 

The mode-coupling piece includes corrections coming from modes other than k in the initial power spectrum. In 
particular, we see from Eq. (|73|) that the one-loop mode-coupling correction involves the convolution of two initial 
power spectra with a specific kernel. At higher loop orders, the mode-coupling pieces involve convolutions of more 
than two initial powers. 

The upshot is that for growing mode initial conditions, we can write the full power spectrum as 

P(k, 77) = G(k, r?)T(fc)T T (fc)G T (k, V )V (k) + P MC (k, 77) • (74) 

At this point we may extend our diagrammatic description: assuming Gaussian initial conditions, the Feynman 
rules for computing the n-loop correction to the power spectrum, SPj? (k, 77) are: 

1. Draw the sum of all tree diagrams corresponding to vp^^k, 77), and the sum of diagrams corresponding to 
v]H 2n- j) (k', 77) as described in §111 Bl including the arrows on the lines. 

2. Formally "multiply" these two sums using distributivity. It is convenient to draw one set of diagrams flipped 
around the vertical axis, so that after the pairing of diagrams for vpw)(k, 77) with those for VE^ 2 ™^ (k, 77), their 
initial conditions face each other. Next, pair the 2n + 2 initial fields into n + 1 initial power spectra (some 
pairings may lead to the same loop diagram). Disregard diagrams that contain a tadpole as a subdiagram. 
(This includes disconnected diagrams in particular.) 

3. Follow Step 131 of §IVCl 

4. Each diagram carries the appropriate factor of 2 r+s inherited form the two original tree diagrams. As noted 
above, it is possible that a specific loop diagram arises more than once. Such a diagram needs to be computed 
only once, but its multiplicity must be taken into account. 

5. The value of 6P£\k,7)) is given by the sum over the the values of these diagrams for all j = 0, . . . , 2n. 

Fig. [5] shows all non-zero diagrams for the power spectrum up to one-loop order, including the symmetry factor 
computed in Step U above. In particular, notice how Eq. (|69[) is very simple to write down by inspecting the second 
and third graphs. The last diagram corresponds to the mode-coupling piece. Evidently Eq. (|73|) may be written down 
immediately by inspecting this graph. 

V. ONE-LOOP CORRECTIONS TO THE POWER SPECTRUM 
A. The non-linear propagator 

In §TVl we defined the non- linear propagator G Q b(k, 77) as the ensemble average of the functional derivative of 
the full solution \I/ a (k, 77) by the initial condition <^°^(k'), without the delta function for momentum conservation 
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(c.f. Eq. I|58p). It was shown to have a perturbative expansion (c.f. Eqs (|59[) and (|60)) ') and terms up to one-loop level 
were computed in Appendix iBl Here we further develop the one-loop expression for the case of N = 2 fluids. 

To begin, let us reconsider Eq. (|62p . we may expand all of the linear propagators, g a b{v): using Eq. (|27p to obtain, 

SG^i^r,) = 4 /"dr/^ ffaMl e Zl ("-"') f d 3 k l7 ^(k 1; k - k,) ^Sccr,^ f dr/' 5> M3 e^'^"> 
Jo f J , Jo , 

t^^)T.9ff'M^ V ''Y.9 g a'M^''V%{k l ) . (75) 



h 

On rearranging the above expression we find 

en 



■ e {li+h-h)r]" 



5G^,(k,r}) = 4 V c' 1 " f dr,'c ih+h - h ^' f dr?"e ( 

x <?a Ml y d 3 ki7^(ki,k-ki) 5cc ^ 2ffde ,^ (76) 

Considering the integrals over 7/ and 77", these may be conveniently evaluated using the algebraic function I2Q2) h, v) 
(c.f. Eq. (|50p V Thus our final expression for the one- loop propagator takes the form: 

«3£>,(k,Tj) = 4 Y, e lr) Hh + h-h,h + h-h;r,) 

h,—,h 

X ffo6.ii/ d 3 kl7bfd( k li k -kl)9cc'J 2 gde,l 3 7efg(-ku (77) 



This expression may be reorganized so that the time dependent algebraic functions, the constant coefficient matri- 
ces g a b,i and the ft-dependent matrices are grouped together. Further, if we take large-scale growing mode initial 

conditions, u a — vffl = (1, 1, 1, 1), then we find that the propagator may be written, 



5G%>,(k,Ti) = 4 ^Hl2 + h-hM + h-h\ri)g a bM 



gcc',l 2 9de,l 3 gff'M gga'.h 



h,-,h 



x Id^^^k^k-kO^^-k^k)^^)^^)^^!) • (78) 

In Appendix ID II we show that the integral over the product of vertex matrices and the power spectrum matrix leads 
to five linearly independent functions of k for each independent product Ti(ki)Tj(ki). Hence, it is possible to write 
the propagator, 



n=l 



where the summation in n runs over the set of independent integrals /* s ° (k) that are defined in Appendix ID II The 
complete results for the Af b s ^(r]), A^ bn (r]) and A^ b s ^(rj) functions are too lengthy to reproduce here, but their general 
structure is the following: 

<bi(v) = E P &iMe ln + Qab.nMV , (80) 
I 

where the summation in I (in this equation only) runs over the set {3, 2, 3/2, 1, 0, —1/2, —1, —3/2, —2, —5/2}, and the 
Pab n i( w 2) an d Qabn( w 2) functions multiplying each 77-structure are at most third order polynomials in W2 (we used 
W\ = 1 - w 2 )- 
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B. Mode-coupling power spectrum 

Consider next the mode-coupling contribution to the power spectrum. At the one- loop level, it is given by Eq. (|73j) . 
We may develop this expression further by again expanding all of the linear propagators using Eq. (|27[) to obtain, 

W^mcM - 2 f dr/ /"V / d 3 k! ^ gac , il e^ f '- f '')7^ e (k 1 ,k-k 1 )^. 9d(iV2 e^X^'^3e Z3 ' ? ' 
Jo Jo J , j i 

x E^/A^-^gj-ki,^ - k) Y,999>M^" E^^ 6 ^">^(fci)^(|k -k!|) . (81) 

On reorganizing this equation and collecting together all of the time dependent terms, we find 
Wj^MoCM) = 2 f dr/ V e ' 1?? e^ +Z3 -' 1 >"' fdrj" V e u "e^+< 6 -^"" 

JO , , 1 JO 111 

x / d^p^^^^k-ki^dd^^ (82) 

The integrals over 77' and 77" are easily performed, and the result can be conveniently written in terms of the I\(l,rf) 
function introduced in Eq. (j4"8)) . Hence, 



5P$uc( k >v) = 2 J2 e hr >h(h + k-h;v) E e u "7i(/ 5 +/ 6 -/4;^) 

Zl^2^3 ^4)^5)^6 

x y d 3 k x ffac,i 1 7ode( k i' k - k i)fdd'^2fee' M^fMlflh (-ki,ki - k)^/ th g h h>,h'P#l> (ki^e-l (l k - k i I ) > ( 83 ) 

which is our final expression for the mode-coupling contribution to the one-loop power spectrum. Upon setting 
large-scale growing mode initial conditions, u a = UeP = (1, 1, 1, 1), the above expression may be written as, 

^mc(M) = 2 y d'kjPoCWdk-kil) 

x E e'^Mfa + is - /u^oci^Cki.k - ki)fl<WM a Se e ',isT ( j/(Ai)T e ,(|k - k x |) 

E e'^Iifo + *6 ~ ^4;r?)56/,u7} 8 i(- k i' k i - K)9 99 >M9hh>MT g '{ki)T hl {\k -k x |). (84) 



x 

^4 ,^5 1^6 



In Appendix ID 21 we show that the integration over the vertex matrices and power spectrum matrices leads to four 
linearly independent functions of k for each independent product Ti(ki)Tj(ki)Tk(\k — ki )Tj(|k — ki |) if (ij) = (kl), and 
six independent functions otherwise. Furthermore, if (ij) ^ (kl), then merely exchanging (ij) ^ (kl) does not lead 
to new independent integrals. Thus, the mode-coupling correction to the one-loop power spectrum can be written: 



sp%mc (k, n) = E (k)A s Z s ° s ° M + f^ sb (k)AZ s ' sb (v) + f* S5S mZ s s W] 

+ D/f W5b (^;f {i w + fr 5b5b (k)Aizf* b ( V ) + f?; sbsbsb (k)Ai^ b ( V )} , ( 85 ) 

n=l 

where the independent integrals fn 638 * 3 ' (k) are defined in Appendix ID 21 The full results for the A^' sks ' (rf) coeffi- 
cients are too lengthy to reproduce here, but their general structure is the following: 

<^%)=£i£#*(«*)«\ (86) 
1 

with I e {4,3,5/2,2,3/1,1,1/2,0,-1/2,-1,-3/2,-2,-5/2,-3} (in this equation only), and Pj[ b %f S ' O2) an (at 
most) fourth order polynomial in W2 (we used w\ = 1 — w%). 
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k [h/Mpc] k [/i/Mpc] 

FIG. 9: Evolution of the initial density perturbations in the CDM plus baryon fluid under the action of the non-linear propagator 
up to one-loop level, as a function of spatial wavenumber k. Left and right panels present results for the CDM and baryon 
fluctuations, respectively. The linear evolution of the initial perturbation in each species has been scaled out, hence all lines 
approach unity on large scales. Solid through to dashed lines show results for redshifts z = {0, 1.0, 3.0, 5.0, 10.0, 20.0}. 

VI. RESULTS: NON-LINEAR EVOLUTION OF CDM AND BARYON FLUID 

A. Propagator 

In Fig.[9]we present the evolution of the initial CDM (left panel) and baryon (right panel) density fluctuations under 
the action of the non-linear propagator, up to the one-loop level in the perturbation series. To show the effects of the 
non-linearity on the dynamics, we have divided through by the linearly evolved initial fields. We see that, for both 
CDM and baryons, on very large scales (k < 0.01 /iMpc -1 ). the modes retain memory of their initial configurations. 
However, on smaller scales the propagator drops away from unity as memory is rapidly lost. The non-linear decay 
scale of the propagator can be represented as the point where G-^ 4>^ / gib^^ — e -1 / 2 with i e {1,3} for CDM and 
baryon density perturbations. We notice that the decay scale becomes smaller as redshift increases, and that for both 
fluid species is essentially the same for all epochs. We also notice that the propagator eventually crosses zero, this is 
unphysical behaviour and signifies the breakdown of the validity of the one-loop expression. 

B. Power spectra of CDM and baryons 

Fig-Hni presents the results for the linear and non-linear evolution of the power spectra in the 2-component fluid 
model. The left and right panels show results for the CDM and the baryons, respectively. As expected from our 
discussion in i jllFI and Fig. [21 we see that the initial CDM power spectrum (zi = 100) is very smooth, whereas that 
for the baryon distribution is highly oscillatory (green solid lines labeled Vc? and T 5 ^? in the plot). As the two fluids 
evolve under gravity the power grows, and in linear theory (thin red lines) larger amplitude BAO are induced in the 
CDM spectrum, whilst those in the baryon distribution are slowly damped. 

The figure also shows how these results change when the non-linear evolution of the density fields is turned on. 
The non-linear power spectra are represented as thick blue lines, and we have included all terms up to one-loop level 
(c.f. Eq. (|7^)1 ). At z — and on very large scales (k < 0.02 h Mpc" 1 ), the linear and non-linear calculations appear 
to agree well (solid lines). However, on smaller scales the agreement begins to breakdown. To investigate this in 
more detail, we take the ratio of the two spectra, and this is plotted in the lower panels of Fig.O From this we 
clearly see that there is a small suppression of power on scales of the order k ~ 0.07/iMpc -1 , this is known as the 
pre-virialization feature (48|. This is followed by a strong non- linear amplification at smaller scales. Interestingly, 
the relative non- linear boost at z — appears to possess almost identical structure for both CDM and baryons. If 
we consider earlier times z > (in the plots denoted by increasingly dashed line styles), then, besides the non-linear 
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FIG. 10: Comparison of the linear and one-loop non-linear power spectrum for the CDM and baryons. Top panels: the 
absolute power as function of wavenumber. Left panels, show results for CDM, and right for the baryons. Thin (red) lines 
represent linear theory, and thick (blue) ones non-linear theory. The solid through to dashed lines show results for redshifts 
z — {0, 1.0, 3.0, 5.0, 10.0, 20.0}. Lowest solid green line shows the initial power spectrum. Bottom panels: show the ratio of the 
one-loop non-linear to linear theory power spectrum. 

boosts being reduced, we note that the curves still maintain the same relative structure. As we will show in i]VII[ this 
can be attributed to the fact that the structure of the leading in 77, one-loop corrections, take the same form for both 
CDM and baryons. 

C. Comparison of the exact 2-component fluid with the approximate 1-component fluid 

We now explore how our exact 2-component fluid results differ from those that would have been obtained had 
we used the 1-component fluid approach, along with the approximate recipe as was described in €JTJ In order to 
concentrate solely on the approximations introduced by using the 1-component fluid approach, we slightly modify the 
procedure for obtaining the initial power spectrum, from that set out in |JJ 

1'. Fix the cosmological model, specifying Sl c and fib, and hence f h . Solve for the evolution of all perturbed species 
using the linearized coupled Einstein-Boltzmann equations. Obtain transfer functions for the CDM and baryons 
at some large redshift say z = 100. 

2'. Use these transfer functions to generate the linear power spectrum matrix in the 2-component fluid theory, at 
the present day: i.e. P a b{k, z — 0) ~ g a a'(z — 0)T a '(k)gbb'(z = 0)Ti,i(k)Ak n . Using this power spectrum matrix, 
define the present day matter power spectrum as 

Pjj(k, z = 0) = (1 - f b ) 2 P s . S c (k, z = 0) + 2(1 - / b )/ b P 5 c 5b (k, z = 0) + (/ b ) 2 /W (k 2 = 0). (87) 

One could further make the approximation that at the present day 8 = (1 — fb)3 c + fbS h ~ hence approxi- 
mating the matter power spectrum as 



j^ A (k,z = o)«iw(k,* = o) 



(88) 
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There rest of the steps are unchanged from those discussed in SjTJ We shall refer to the initial conditions obtained using 
Eq. (|87[) as the "exact", and those produced using Eq. ([55]) as the "approximate" 1-component fluid initial conditions, 
respectively. 

Figs[TT] and rr2J present the ratio of the exact linear and non-linear CDM (Fig. [TTj) and baryon fFig.[T2"|) power 
spectra, with the power spectrum obtained from modeling CDM and baryons as an effective single matter fluid. The 
left and right panels show results obtained using the exact and approximate initial conditions for the 1-componcnt 
fluid calculation. In all cases, the non-linear evolution is followed using the the RPT framework, including all terms 
up to the one-loop level. For the case of CDM, we clearly see that at high redshift the approximation is rather poor, 
there being significantly more amplitude in the exact N = 2 model than the effective N = 1 model. We also note 
that there are strong BAO features in the ratio, and these increase with increasing redshift. These features are due 
to the fact that the z = transfer function is used to model the initial conditions of the N = 1 matter fluid, and 
therefore can suppress power and artificially enhance the BAO. At early stages in the evolution, the systematic errors 
are between ~ 4 — 5% at z = 20, and between ~ 2 — 3% at z = 10. However, at later times the approximation is 
much better, and by z — 3 the errors are < 1%, and by z — are < 0.5%. 

On the other hand, the case for the baryons is worse. We see that the exact N = 2 results show less power than the 
approximate N = 1 counterpart and that this suppression increases with increasing redshift. Again we notice that 
there are BAO features in the residual, and the oscillation amplitude is much more significant. At early stages in the 
evolution the systematic errors are large: at z — 20 they are ~ 25%, and by z = 10 are ~ 15%. At later times the 
approximation is still quite poor: at z = 3 the errors are of the order ~ 5%, and by z = are still between ~ 2 — 3%. 
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FIG. 11: Ratio of the CDM power spectrum obtained from the CDM and baryon 2-component fluid model, to the total matter 
spectrum obtained from the 1-component fluid approximation, as a function of wavenumber. Thin (red) lines show the ratio 
of power spectra in linear theory, while thick (blue) lines show the ratio of one-loop non-linear power spectra. Left and right 
panels: show results for the exact and approximate 1-component fluid initial conditions, respectively. Again, the solid through 
to dashed lines show results for redshifts z = {0, 1.0, 3.0, 5.0, 10.0, 20.0}. 
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FIG. 12: Ratio of the baryon power spectrum obtained from the CDM and baryon 2-component fluid model, to the total matter 
spectrum obtained from the 1-component fluid approximation, as a function of wavenumber. Thin (red) lines show the ratio 
of power spectra in linear theory, while thick (blue) lines show the ratio of one-loop non-linear power spectra. Left and right 
panels: show results for the exact and approximate 1-component fluid initial conditions, respectively. Again, the solid through 
to dashed lines show results for redshifts z = {0, 1.0, 3.0, 5.0, 10.0, 20.0}. 



23 



1.0003 



1.0002 



1.0001 




0.001 0.01 0.1 

k [h/Mpc] 



1.0 



0.998 



0.996 



0.994 







_ 


\ 




\ / 1 1 I 19 / i/y^-^T ™ 


- 


V/ V ^^ N 








V \ . 











0.001 



0.01 



0.1 



k [/i/Mpc] 



FIG. 13: Ratio of the total matter power spectrum obtained from the 2-component fluid model for CDM plus baryons, to 
that obtained from the 1-component fluid approximation, as a function of wavenumber. Again, thin (red) lines show the ratio 
of power spectra in linear theory, and thick (blue) lines the ratio of one-loop non-linear power spectra. The solid through to 
dashed lines show results for redshifts z = {0, 1.0, 3.0, 5.0, 10.0, 20.0}. See end of ijVIDI for an explanation as to why the linear 
theory ratios are independent of time 



D. Impact on total matter power spectrum 



Since some cosmological probes are only sensitive to the total mass distribution, such as weak gravitational lensing, 
it is interesting to ask, what is the error between the exact N = 2 and the approximate N = 1 fluid dynamics, incurred 
when modeling the total matter power spectrum. Recall that the total matter power spectrum is defined as, 

P rs (k, z) = (l- f h ) 2 P g ^ (k, z) + 2(1 - f h )f h P s , sb (k, z) + (f h ) 2 P 5bsb (k, z) . (89) 

The results obtained by using the exact and approximate initial conditions for the 1-component fluid model are 
plotted in the left and right panels of Fig. Q]|] respectively. We see that the single-fluid approximation based on exact 
initial conditions is excellent, being accurate to < 0.03% for all times and scales considered. As we will show in Will 
this spectacular success can be attributed to the fact that the leading in 77, one-loop corrections, depend only on the 
total mass distribution. 

Using the approximate initial conditions for the single-fluid calculation still leads to very good agreement with the 
full 2-component fluid result on large scales k ~ 0.2 /iMpc -1 , with < 0.5% accuracy for all the times considered. On 
smaller scales the approximation becomes worse, being of order ~ 0.7% at k ~ 1.0/iMpc . However, this departure 
is unlikely to be accurate, since we can not trust the PT results on these small scales. This owes to the fact that 
the missing higher order loop corrections become increasingly more important. Nevertheless, this point is rather 
academic, since using the exact initial conditions will give better than 1% precision for the non-linear matter power 
spectrum down to scales of the order k ~ 1.0 /iMpc -1 , as required by future weak lensing surveys, such as Euclid. 

Finally, we note that one can easily show that the linear mass power spectrum simply scales with the square of 
the linear growth factor, even in the 2-componcnt fluid theory. Hence, the ratio of the 2-component fluid linear mass 
power spectrum to the 1-component fluid model is constant in time. Furthermore, this constant is simply unity, if 
the 1-component fluid model is set up with exact initial conditions. 



VII. APPROXIMATE DESCRIPTION OF THE MULTI-FLUID DYNAMICS 



Given that the full expressions for the one-loop corrections to the power spectrum are quite formidable and time 
consuming to compute, it is of use to attempt to develop an approximate description of these expressions. We present 
this below. Our approximation scheme is based on insights concerning the late-time behaviour of the one-loop power 
spectrum, and it is valid for growing mode initial conditions and, although here we present the case of a 2-componcnt 
fluid, can be generalized for A^-component matter fluids. 
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A. Approximate form of the reducible correction 

To begin, we recall that the one-loop reducible correction to the power spectrum reads (c.f. Eq. (|69|) ) 

wi&M = Lc(r/)5Gi 1 d ) (k,r ? )+5GW(k,r ? ) 56d (r ? )l P { $(k). 



(90) 



In the limit of late times, rj — > oo, we can discard all but the fastest growing modes of the linear propagator, g a b, and 
one-loop correction to the propagator, SG^ b . For the linear propagator we have simply gabiv) ~ * 9ab.i eTI , while for 
the one-loop correction we find 

,5G$(M)- 



3 Wl f(k) 2 Wl f(k) 3w 2 f(k) 2w 2 f(k) 

3wig(k) 2w\g{k) 3u>2<?(fc) 2w2g(k) 

3wxf(k) 2 Wl f{k) 3w 2 f{k) 2w 2 f(k) 

3w\g(k) 2w\g{k) 3w2g(k) 2w2g(k) 



(91) 



where we have 
f(k) = 

g(k) = 



d 3 k! 

504fc 3 fcj i 

d 3 k! 

168fc 3 fcf 



6k 7 h - 79fc 5 fc 3 + 50k 3 kf - 2\kk\ + - (fc 2 - fc 2 ) 3 (2fc 2 + Ik 2 ) In j* ^ 

4 |fc + fci| 2 

6k 7 h - 41fc 5 fc 3 + 2k 5 k\ - 3kk 7 + -(fc 2 - fc 2 ) 3 (2fc 2 + fc 2 ) In lf~M' 
1 1 1 4 V 17 1 |fc + fci| 2 



( 92 ) 



(93) 



These are just the / and g functions of (27j, and vW{k) is the initial power spectrum of the total matter fluctuation, 
8{k) = wi8x(k) + w 2 5 2 (k): 

r^(k)^wjv[ 1 \k) + 2w 1 W2V[ 2 ) (k)+w 2 2P^ > (k) = [w 1 T 1 (k) + w 2 T 2 (k)] 2 V (k), (94) 

where we have assumed growing mode initial conditions. Upon evaluating Eq. (|90[) with the approximate form of the 
propagators as given above, we obtain 



2/(fc) f(k)+g(k) 2/(fc) f(k)+g(k) 

f(k)+g(k) 2g(k) f(k)+g(k) 2g(k) 

2/(fc) f(k)+g(k) 2/(fc) f(k)+g(k) 

f(k)+g(k) 2g(k) f(k)+g(k) 2g(k) 



For early times on the other hand, as rj — > we have simply that 6P^ R (k, rj) — > 0. In fact, explicit computation 



(95) 



shows that in this limit, the reducible correction vanishes as S P^ R (k, rj) ^oc (e 1 ' — l) 2 . Since the vanishing of 



SP^ b R (k, ?y) as rj — > is required on physical grounds, we would like to implement it even in the approximate 
description. To do so, we will simply change the time dependence of Eq. ((95)) as e iri — > e 2r> (e n — l) 2 . Thus we define 
the approximate one-loop reducible correction to the power spectrum as 



< ) R,A(k,')) 



2/(fc) f(k)+g(k) 2/(fc) f(k)+g(k) 

f(k)+g(k) 2g(k) f(k)+g(k) 2g{k) 

2/(fc) f(k)+g(k) 2/(fc) f(k)+g(k) 

f(k)+g(k) 2g(k) f(k)+g(k) 2g{k) 



vfMk)^-if 



(96) 



Needless to say, the fastest growing mode is not affected by this substitution. 



B. Approximate form of the mode-coupling correction 

Keeping only the fastest growing mode of the one-loop mode-coupling correction to the power spectrum, we find 



6 P. 



(i) 



afc,MC 



(M) 



2F(k) H(k) 2F(k) H{k) 

H(k) 2G(fc) H(k) 2G(k) 

2F(k) H(k) 2F{k) H(k) 

H(k) 2G(fc) H(k) 2G(k) 



(97) 
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FIG. 14: Ratio of the approximate to exact one-loop non-linear matter power spectra for CDM (left panel) and baryons (right 
panel) in the 2-component fluid model for CDM plus baryons. As before, the solid through to dashed lines show results for 
redshifts z = {0, 1.0, 3.0, 5.0, 10.0, 20.0}. Approximations are good to better than 1% up to k « 0.2/iMpc -1 . 



where we have defined 

fc 4 [7kx + h (3 - 10.x ; 



G(fc) 
H(k) 



d^k. 



d 3 ki 



196fc 2 (fc 2 - 2fcfcia; + fc : 
fc 4 [7fcx-fci(l + 6a; 2 )] 



_ v (o) 



(h)V 



(0) 
55 



d 3 ki 



196fc 2 (fc 2 
fc 4 [(60x 4 



2kk\x 
8x 2 - 



v (o) (k )v (o) 

2 ' 55 ^ X '' 55 



fc 2 - 2fcfci.T + k\ 



fc 2 - 2fcfcia; + k\ 



3) fc 2 - Ukk lX (8x 2 - 1) + 49fcV] (0) (0) 

' 55 V Kl ^55 



98fc 2 (k 2 ~ 2kk x x + k\y 



2kk\x 



kj 



(98) 
(99) 
(100) 



where x is the cosine of the angle between k and ki, i.e. k • ki = kk\x. 

For early times, the situation is completely analogous to the case of the reducible correction: as r\ 
that 5P}& Mc(k, rj) —> and in fact, S P^ MC (]s., r]) ^oc (e v — l) 2 in this limit. So again we implement 
of the mode-coupling correction as rj — > by simply changing the time dependence of Eq. (|97|) as e 4l? — 
Then the approximate one-loop mode-coupling correction to the power spectrum is defined as 



SP. 



(i) 



ai,MC,A 



CM) 



2F(k) H(k) 2F(k) H(k) 

H(k) 2G(k) H(k) 2G(k) 

2F(k) H(k) 2F(k) H(k) 

H(k) 2G(k) H(k) 2G(k) 



(e" - 1)' 



— > we have 
the vanishing 
-e 2 "(e" - l) 2 . 



(101) 



Clearly the fastest growing mode is not affected by this substitution. 



C. Comparison of the approximate and exact 2-component fluid models 

Fig.HU shows the ratio of the approximate non-linear CDM (left panel) and baryon (right panel) power spectra 
with the exact 2-component fluid non-linear power spectra: 

<a(M) _ iff (fc, z) + 6P$ R A + SP^ MC A (k, z) 
P%\k,z) ' P! l l ) (k,z) + 6P$ R + 6P$ MC (k,z) ' 

where the first term in the numerator is given by Eq. (jfj?)) with g a b being the full linear propagator. For the case of 
CDM, the agreement is very good, being < 0.2% for fc < 0.2/iMpc" 1 (i.e. on scales large enough for one-loop PT to 
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be trusted), and for all times considered. The approximation performs worse for baryons, however, the agreement is 
still quite good, being < 0.7%, on scales of up to k ~ 0.2 /iMpc" 1 , again for all considered times. 

We note in passing that the success of the approximation - which was based on the late time behaviour of the 
one-loop correction - for relatively early times (say z = 10-20) is simply a consequence of the fact that at early times 
the one-loop correction to the power spectrum is small. Hence, even if the approximation of the one-loop correction 
itself is not very accurate at such times, this does not influence the fact that the full non-linear power spectrum is 
still well approximated. Clearly, if all of the SP^s are small compared to P(°) in Eq. (|102p . then the ratio is close to 
one, irrespective of whether the individual approximations <SPW rj dP»P are particularly good or not. 

We thus conclude that the one-loop corrections to the power spectrum for both CDM and baryons are very well- 
modeled by the leading late time terms, for all scales and times considered. Incidentally, this observation explains two 
phenomena we encountered earlier. Firstly, as noted in i jVI Bl the structure of the non-linear to linear power spectra 
for the CDM and baryons (see bottom panels of Fig. [T0|) are almost identical. This is because the leading late time 
corrections in Eqs (|96[) and (jlOip for both CDM and baryons are seen to take precisely the same form. Secondly, the 
spectacular success of the single-fluid approximation in modeling the total mass power spectrum, noted in WIDI (see 
Fig. [T3")) . is due to the fact the leading late time terms depend only on the average initial field, and not on the two 
fields separately. Thus, treating the total 2-component fluid as an effective 1-component fluid is exact at the level of 
the leading non-linear corrections. 

Finally, this approximate theory is of great practical advantage, since it reduces the number of numerical integrals 
to be evaluated by roughly a factor of 10 for the case of the 2-component fluid, so giving an order of magnitude speed 
up in computational time. Further, the approximate theory can be represented by a compact set of equations. 

VIII. SUMMARY AND CONCLUSIONS 

In this paper we have generalized the matrix based RPT formalism of Crocce and Scoccimarro [l| to the problem 
of dealing with TV-component fluids, each of which evolves from a distinct set of initial conditions and contributes to 
the matter-density and velocity fields. 

In ||TT]we showed that, for TV-component fluids, the essential structure of the RPT framework remained the same 
as for the TV = 1 case - that is, the equations of motion could be solved in exactly the same way. However, the details 
of the building blocks of the theory did change: i.e. the vertex, propagator and initial conditions. For TV = 2 the 
vertex matrix offered no surprises. The linear propagator was more interesting, and we showed that besides the usual 

growing and decaying modes {-D+ oc e'', D_ oc e~ 3 ''/ 2 }, there were two additional eigenvalues: a static and decaying 

mode {Astat oc const, D_ oc e - ''/ 2 }. Interestingly, in going to TV > 2 no further time dependence was gained. Unlike 
the 1-component fluid case, it was no longer possible to choose pure growing modes, unless all of the fluids evolved 
from identical initial conditions. However, we showed that one can still set up the initial conditions to possess pure 
large-scale growing modes. All our results were obtained for this specific choice of initial conditions. 

In mill it was shown that the equations of motion for the TV-component fluid case could be solved using a power 
series solution and that the solutions could be recovered order by order. A graphical representation of the solutions 
was also described. 

In mV\ we discussed the statistics of the fields: in particular the covariance matrix of the TV-component fluid 
perturbations, and the fully non-linear propagator. This latter quantity can be understood to be the statistic that 
provides information on the memory of the evolved field to the initial conditions. It also possesses a perturbative 
expansion and we gave expressions up to one-loop order. A diagrammatic representation was also discussed. 

As for the case of the 1-component fluid RPT, we then showed that for the TV-component fluid RPT the full non- 
linear power spectrum also possessed a series expansion, and that the series could be grouped into two sub-series, 
"reducible" and "mode-coupling" . It was shown that the reducible terms were proportional to the initial power 
spectrum and non- linear propagators. Again we gave complete details for the series up to the one- loop level in the 
expansions. A diagrammatic description was also discussed. 

In |JV]we developed further the one-loop expressions for the propagator and the mode-coupling contribution to the 
power spectrum. These expressions were then given for the explicit case of a 2-component fluid. 

In § VII we applied the formalism to the problem of modeling the non-linear evolution of a CDM and baryon fluid 
in the ACDM paradigm. This enabled us to test the validity of the approximation of treating CDM and baryons as 
an effective 1-component fluid evolving from a single set of initial conditions, as is currently standard practice. 

For the case of CDM, it was shown that the approximation was very good for z < 3, with the exact and approximate 
CDM power spectra differing by < 1%. However, for higher redshifts the approximation became progressively worse, 
it being of the order ~ 3% at z = 10, with the power in the exact theory being amplified and having weaker BAO 
features than the approximate theory. 
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For the case of baryons the situation was worse: at z = 20 the exact and approximate spectra differed by ~ 25%, 
and at z = 1 by ~ 15%. At later times the approximation was still quite poor, and at z = 3 the errors were of the 
order ~ 5%, and by z = were still between ~ 2 — 3%. The power in the exact theory was suppressed in amplitude 
but possessed stronger BAO than the corresponding power for the approximate theory. 

These conclusions remained essentially unaffected, when the so-called approximate 1-componcnt fluid initial condi- 
tions were employed for the effective 1-component fluid computation. 

Lastly, we computed the total non-linear matter power spectrum and found that the exact and approximate theory 
agreed to within < 0.03%, when the so-called exact 1-component fluid initial conditions were used. Employing the 
approximate 1-component fluid initial conditions on the other hand leads to an agreement to within < 0.5% on scales 
where the perturbation theory was valid k < 0.2/iMpc -1 and for all times of interest. Deviations on the order of 
~ 0.7% were found on smaller scales with this approximation. 

Our main conclusions therefore are: 

• For theoretical modeling of the low-rcdshift CDM distribution, the approximate 1-component fluid treatment 
provides a good approximation. For higher redshift studies, such as those that probe the epoch of reionization, 
or attempt to model the first objects that form (stars/haloes), then one must be careful to use the appropriate 
CDM transfer function for the redshift of interest. Otherwise, systematic errors will be present at the level of 
several percent. 

• For theoretical modeling of baryons in the Universe, a 1-component fluid treatment leads to significant > 1% 
errors at all times. This implies that cosmological probes that are primarily sensitive to the distributions of 
baryons in the universe, such as: 21 cm radiation from neutral hydrogen, or the Lyman alpha forest absorption 
lines in high redshift quasars, can not be modeled using dark matter only simulations at high precision. Instead 
the baryons must be modeled using a 2-component fluid system of CDM and baryons, with the initial distribu- 
tions being specified by the linear Einstein-Boltzmann solutions. The correct treatment of CDM and baryon 
initial conditions may also affect how galaxies form. 

• Observational probes for cosmology that are sensitive to the total mass distribution, such as weak gravitational 
lensing, can be accurately interpreted using an effective single CDM plus baryon fluid. Thus current modeling 
technology is good enough for high precision work, but it is highly desirable to use the exact prescription when 
specifying the initial conditions for the effective single fluid. 
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APPENDIX A: A RELATION INVOLVING THE CONFORMAL EXPANSION RATE 



We prove that dH/dr = —H 2 [fi m (a)/2 — Qa(g)]. To begin, we recall that for a homogeneous and isotropic universe, 
the Einstein field equations reduce to: 



^ = -4*G(3p + p); a 2 + K= 8 -^ 
a 3 



(Af) 



where dots denote differentiation with respect to cosmological time t. Then for the ACDM model with matter density 
P = Pm + PA and pressure p = p\ = —p\, we find 



a 
a 



4nGp m 8irGpfi 



3 



= -H 2 (a) 



fim(a) 



Consider now, 



dr dr 

The last two equations together give the desired result 

dr 



dH d dt d . . . 



dr dt 



n m (o) 



(A2) 



(A3) 



(A4) 



APPENDIX B: CALCULATING THE PROPAGATOR 

1. The non-linear propagator 

In this appendix we provide explicit details of the calculation of the propagator for the PT solutions up to third 
order in the expansion for ^i™" 1 . To begin, recall that the propagator is defined through Eq. (|55)) . as 



4 0) (k') 



(Bl) 



G ab (k, v )5 D (k-k') EE 

If we insert the series expansion for ^i"^ as given by Eq. (|41|). then we obtain a series expansion for the propagator 



G ab (k,r 1 )6 D (k-k') = Y, 



, ^i 0) (k') 



(B2) 



The next step is to calculate the functional derivatives of the perturbative solutions. Functionals are essentially 
functions that depend on the full shape of another function over a certain domain. We shall denote functionals using 
the notation, F[y(x)], and this has the specific meaning, 



(B3) 



F [y(x)} = / dxg(y(x)) . 

J a 

where g(y) is some arbitrary continuous function. Functional derivatives may be taken using the relation, 



dF 



da; 



8F 



5y(x) 



8y{x) 



(B4) 



y°{x) 



where in the above we identify 8F / 5y{x)\ y0 ^ as the functional derivative about the function, y°(x). 



2. Propagator for solutions: n — {0, 1,2} 
To see how we may use Eq. (|B4[) . let us explicitly apply it to the solutions <j("), for n = {0, 1, 2}. 
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• Solution for n = 0: 

For the linear theory solution we have, ^i ' (k, 77) = <7a6 ( 7 7)0t,°' ) (k) . This is not in the form of Eq. (|B3j) . but we can 
easily write it in such a form using the Dirac delta function, 

*i°)(k,77)[0i O) (k')] = J dW(k - k')g ab ( V )^\k') . (B5) 
Now consider a small change Stfif 1 ^ (k') to the initial conditions 4>jP\k'), whereupon we have 
*<°>(k,»7)[^ 0) (k') + <S^ 0) (k')] = / dW(k - k')g ab { V ) \<j>f\k') + <5^ 0) (k') j 

(B6) 



= *(°)(k,r,)[0i O) (k')]+ Jd 3 k>8 D (k-k>)g ab (r))5$\k'). 
Taking the first term on the right-hand-side over to the left, we find 

d*i°)(k,, 7 )[0i O) (k')] = /dW(k-k'). 9Q6 (r/)^i 0) (k') . 



and on comparing the above with our definition from Eq. (|B4[) . we identify the functional derivative as, 

<s*£ 0) (M) _ sDt 



WW 



= S D (k-k')g ab (r 1 ) . 



(B7) 



(B8) 



Finally, we are interested in the ensemble average evolution of the mode, and since there are no initial fields remaining 
there is no stochasticity left in the relation and we simply have, 



= S n (k-k')g ab (r,) . 



(B9) 



• Solution for n = 1: 



Let us now turn to the more complicated situation of computing the functional derivatives of the first order solution. 
The first order solution is a functional of the form, 



*( 1 )(k, 7? )[#(k 1 ),^? ) (k 2 )] = ^ difg^n-tf) J d 3 k 1 d 3 k 27 W(k,k 1 ,k 2 ) 5cc ,(r?O^ ) (ki)w(^)0, ( ° ) 



(k 2 ) . (BIO) 



Repeating the procedure of above, we consider a small change (k) to the initial conditions cjyf, (k), whereupon 
we have 

d*PM[$\</>W] = j\rf gab (n-rf) |d 3 k 1 d 3 k 27 W(k,k 1 ,k 2 ) ffcc ,( ? 7')^(»?') 
X '^(k!)^(k a ) + <$\k x )84>f{k 2 ) 

= 2j V d V ' 9ab (r 1 -r,') J d 3 k 1 d 3 k 27 ^(k,k 1 ,k 2 ).g cc ,(r/)^? ) (k 1 ).g^(ry')^? ) (k 2 ) , (Bll) 

where the last equality follows from the fact that owing to the symmetry of the vertex matrix 7 ^(k, ki,k 2 ), the 
solutions are symmetric to changes {ki,k 2 } — > {k 2 ,ki} and {c, d} — > {d,c}. 

We would like to take the functional derivative with respect to the arbitrary initial condition 5(f>^}(k'), and so we 
use the relation 

6^(k t ) = J d 3 k J -^(k i - k^frf^k,) , (B12) 
to rewrite the initial conditions in terms of the field that we vary. Hence, Eq. (|B11[) , becomes 

d*^ (k, » 7 )[^ ) (kx), 0^ (k 2 )] = J\ V 'g ab ( V - V ') J d 3 k 1 d 3 k 27 ^(k,k 1 ,k 2 ) 3cc ,(r/)^? ) (k 1 ).g^(ry') 

x Jd 3 k'5 D (k 1 -k')Sf, a ,S<f> { °\k') . (B13) 
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On comparing the above expression with Eq. (|B4[) , we identify the functional derivative as 
M X) (M) _ xDl 



d D (k 1 -k') J V 6rfgab(v-rf) J d 3 k 1 d s k 2 7^(k,k lj k 2 ) Scc> (»/)^ 0) (ki)(; do ,(»/) 



(B14) 



On taking the ensemble average we have ((f>£!\ki)J = 0, and hence 

/ M 1) (k,7 7 ) \ 

\ / 



(B15) 



One direct corollary of this result is that, the propagator for the perturbative solutions which involve an even number 
of initial conditions is zero. Hence, only odd powers of <^ ' contribute. 

• Solution for n = 2: 

Based on the previous calculation, we expect that the first non-zero contribution to the non-linear propagator is n = 2, 
since it involves three initial conditions. This can be computed as follows. The solution is a functional of the form, 

*( 2 )(k, J7 )[^? ) (k 1 ),^ ) (k 3 ),^ ) (k 4 )] = 2 fdrfg^n-rf) f d^d^Q^kxM^M^^rf) 



= 2jfWfl«&fa-»fl J d 3 k 1 d 3 k 27 W(k,k 1) k 2 ) 5cc ,(r ? ')^? ) (ki; 
x J d V "g de ( V ' - r,") J d 3 k 3 d 3 k 47 ^(k 2 ,k 3 ,k 4 ) 
x 9ff>{rf')4>f^)9 B Av[')4>fQ*) , 



(B16) 



where the factor of two comes from the fact that the first equation is symmetric to interchanging, {c, d} — > {d, c] 
and {ki, k 2 } — » {k 2 , ki}. Repeating the procedure of above, we now consider a small change Scjyf, (k') to the initial 



conditions 0^, (k'), whereupon 
d*l 2 )(k, ?? )[^ ) (k 1 ),^ ) (k 3 ),^(k 4 )] = 2j n dr 1 'g ab (r 1 -r ] ') J d 3 k 1 d 3 k 27 ^(k, k 4 , k 2 )g cc ,(v') £ d V " g de ( V ' - V ") 

x f d 3 k 3 d 3 k 47 ^(k 2 ,k3,k4). g// ,(ry"). 9 ^( ?? ") [^? ) (k 1 )4° ) (k 3 )^? ) (k 4 ) 



+ ^? ) (k 1 )^(k3)^(k4) + ^(k 1 )^(k 3 )^(k 4 ) 



(B17) 



We notice that the last two terms are symmetric to changes {k 3 ,k 4 } — > {k 4 ,k 3 } and {/,<?} — » {.9,/} and this gives 
us a factor of 2. Next we use Eq. (|B12I) to repeatedly rewrite the perturbations Scf) 1 - ^ in terms of the the arbitrary 

initial field dcfr^) (k'). Hence, 

dvp( 2 )(k,r/) = 2 r dv'g ab (v-v') I d 3 k 1 d 3 k 27 ^(k,k 1 ,k 2 ) 5cc ,(V) f d^'gdeirj'-v") 
Jo J Jo 

x J d 3 k 3 d 3 k 47 ^(k 2 ,k 3 ,k 4 ) 5// ,( ? 7")5 99 '('7") 

>(k! - k')^ o/ 0W(k3)4 O) (k 4 ) + 2^(k 1 )4°)(k 3 )^(k 4 - k')^l <5^(k') . (B18) 



x / d 3 k' 



On comparing the above with our definition for the functional derivative as given by Eq. (IB4|1 . we may make the 
following identification: 



# } (k') 



1 1 -2/ dV5a6(»7-»7')/d 3 k 1 d 3 k 27 W(k,k lj k 2 )^(»/)/ d,V(>?'-l") 
Jo J Jo 

| d 3 k 3 d 3 k 47 5 g (k 2 ,k 3 ,k 4 ) 5/r (^")^( ? 7") 
6 D (k 1 - k^^ik^ik,) + 2^? ) (k 1 )^ ) (k 3 )5 D (k 4 - k')5 9 \, 



(B19) 
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Finally, we take the ensemble average and find that the term in square brackets becomes, 



^( kl - k'fe (^(k 3 )<^(k 4 )\ + 2 (^(k^Ck,)) S °(k 4 _ V )5 K a , 



5^(ki - k')^ a ,P^(k 3 )^(k 3 + k 4 ) + 2^(k 1 )«5"(k 1 + k 3 )^(k 4 - k')e o 



(B20) 



Consider the first term in this bracket, we note that it contains a S D (k 3 + k 4 ). If we integrate this first term over 
k 4 then the resulting expression involves 7 g^ g (k 2 , k 3 , — k 3 ), which from momentum conservation means that waves 
with equal and opposite momenta enter the interaction and so completely cancel and no wave exits. Thus, the above 
expression reduces to 



'<S*i 2) (M) 



= 4 fdrfg^r) - tf) [ d 3 kid 3 k 2 7« (k l5 k 2 )5 D (k k, k 2 )g cc ,(rf) f tof'gtJrf - v") 
y(k') /Jo J Jo 

x J d 3 k 3 d 3 k 4 7^(k 3 ,k 4 )^(k 2 -k 3 ^ . (B21) 



If we integrate over k 4 , then k 3 , and finally over k 2 , then we obtain the desired expression, 

/**< 2) (M)\ _ xD 



WW 



«^(k-k')x4 fdrfgahiri-rf) f d 3 k 1 ^(k ll k-k. 1 )g cc ,{r/) f d V " g de ( V ' - r,") 
Jo J Jo 

x tl(-^)9fAri' , )9 a u'{r 1 ' l )V^} l {k l ) . (B22) 



Finally, if we assume large-scale growing mode initial conditions u a = (1,1,. ..,1,1), then we find the first non-zero 
perturbative correction to the propagator is given by: 



«?£>(M) = 4 r drf g ab (ri-rf) [ d 3 ^^ ( kl , k - ^g^T^kx) f d V " g de (r,' - V ") 
Jo J Jo 



x 7 W (-k x , k')g ff , {rf')T f , (h )g ga , ( V ")V (k x ) 
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APPENDIX C: CALCULATING THE MODE-COUPLING POWER SPECTRUM 

In this appendix, we present details of the computation of the mode-coupling contribution to the one-loop power 
spectrum given by Eq. (|73|) . In Eq. (|68p we identified this term as, 



*W(k)*«(k')) = dr n g ab (v-vi) Id^id^^C^k^^^^O^^o^Ck!)^ ^^) 

x J\ m g a , v { V - m ) J d 3 k 3 d 3 k 47 ^, d ,(k',k 3 ,k 4 )^ e -(r72)5d'/'('72)^? ) (k 3 )^° ) (k4) 

= Jf[ d3 ^J o di] 2 ga>b>(v - mWvdd'^y ^^)9c'e-{m)9d'f'{m) 



i=l 



X / d7 7lffa6 (7 7 - 7 7l ) 7 W(k,k 1 ,k 2 ) ffce (7 7l )^(7 7l ) (4°)(k 1 )^(k 2 )^? ) (k 3 )^(k 4 )) .(CI) 



On assuming Gaussian initial conditions we may use Wick's theorem to rewrite the product of initial fields, 

^(kO^CkaJ^fka)^^)) = [Vef{^i)5 D ^i +k 2 )P e , f ,(k 3 )S D (k 3 + k A ) 

+ V ee/ (k 1 )5 D (k 1 +k 3 )V fr (k 2 )5 D (k 2 +k 4 ) 
+ P e/ ,(k 1 ) < 5 D (k 1 +k 4 )^(k 2 ),5 D (k 2 + k 3 )] . 



(C2) 



We notice that the first term in this structure involves Dirac delta functions with 8 D (ki + k 2 ) and S D (k 3 + k 4 ), if we 
were to perform any of the fc-integrals, this would lead to a vertex matrix 7^<j(k, ki, — ki) = 0. Thus from momentum 
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conservation the first term vanishes. The next two terms are non-zero, but are symmetric to relabelings, and so give 
a factor of 2. Thus we have, 



*? ) (k)*i/ ) (k / )) = 2j dnig a b{n-m)9ce(rii)g <v (ni) J d^g^b-iv-^go'^^gapim) 
r 4 

x / J]d 3 k i7 ^(k,k 1 ,k 2 ) 7 ^ M ,(k',k 3 ,k 4 )P ee ,(k 1 )^(k 1 + k 3 )V fr (k 2 )6 D (k 2 + k 4 ) . (C3) 

i= 1 

Considering the fc-space integrals, if we perform the k4 and IC3 integrals then we arrive at, 

- | d 8 k 1 d% 7 ^(k,k 1 ,k a ) 7 W d/ (k , ,-k 1 ,-k 2 )7' ee .(k 1 )7' // ,(k 2 ) . (C4) 

If we now use 7^ c (k, ki, k 2 ) = 7^; c (ki,k 2 )<5 z:, (k — ki — k 2 ) to factor out the Dirac delta functions from the vertex 
matrices, then we find 



(^(W^i'V)) = 2 y dr] 1 g ah {i 1 -r] 1 )g ce {r}i)g df {r 1 i) g a , b ,{ v - m ) g c , e ,( m ) g d , f ,(r] 2 ) 

x /" d^id^^ki.ka^k-k!-^ (C5) 

On computing the integral over k 2 and factoring our the delta function S D (k + k'), we find that our general expression 
for the one-loop mode-coupling contribution to the power spectrum is given by, (c.f. Eq. ([73")) ) and see also [l|) 

5P aa>, Mc( k ' ? /) = 2 / d3k i / d Vi [ d?7 2 5 ab (?7-?7i) 7 ^(ki,k-ki).g ce (?7i) 5d/ (?7i) 
J Jo Jo 

x g a 'v (r) - r? 2 )7bvd'(- k i 5 k i - k ) 9c e ' (»») gap {r\i)V^ (k^TV (|k - k x |) . (C6) 



Finally, on choosing large-scale growing mode initial conditions, wi = (1, . . . , 1), the above expression reduces to 

€,mc(M) = 2 1 d 3 k 1 7'o(fci)7'o(|k-k 1 |) 

x / di7iSf o6 (rj-7?i) 7 ^(ki,k-k 1 )fli ce (j7i)Te(A;i)sf4f(7yi)T/(|k-ki|) 
Jo 



x / d?7 25a / 6 /(77-?7 2 ) 7 ^, d ,(-ki,ki -k).g c / e ,(?7 2 )T e /(fci) 5d / / /(7] 2 )T / /(|k-ki|) . (C7) 
Jo 

APPENDIX D: SOME MOMENTUM INTEGRALS 
1. Momentum integrals arising in computing <5G^(k, rf) 



When computing 5GQ (k, 77), we encounter the momentum integral 

J d 3 k 1 7^(k 1 ,k-k 1 )7g s (-k 1) k)^,°),(fc 1 ). (Dl) 
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Since all non-zero matrix elements of 7^ c (ki, k^) are either a(ki,k2)/2, a(k.2,ki)/2 or /3(ki,k2), we see that only 
nine distinct types of integrals can arise in the computation. These are 



lfV(k) = y'd 3 k 1 a(k 1 ,k-k 1 )a(-k 1 ,k)P. i 
d 3 ki a(ki, k - ki)a(k, -kx)P,, 
d 3 ki a(k - ki , ki)o ! (-k l! k)V, 
d 3 ki a(k - ki, ki)a(k, -ki)7> 



l( si (k) 



lf si (k) 



lf si (k) 



d 3 k 1 /5(k 1 ,k-k 1 )a(-k 1 ,k)P 1 



l( 5i (k) = J d 3 k 1 a(k 1 ,k-k 1 )/3(-k 1 ,k)P 1 



(0 

ij 

l( 5 \k) = J d^alk-k^kO^-k^k)^ 

l( 5 \k) = I d 3 k 1/ 3(k 1 ,k-k 1 )a(k,-k 1 )^° 

,(o 



(fci) 
(fci) 
(fci) 
(fci) 
(hi) 

(fci) 
(fci) 
(fci) 



1 



-1 d^^^x) 



(0), 



d 3 ki 
d 3 ki 
d 3 ki 



fc 2 + fc 2 (fc 2 -fc 2 ) 2 \k-h\ 



4fc 2 



In- 



\k + h 



(D2) 
(D3) 

^(fci), (D4) 



3fc 2 -fc 2 (fc 2 -fc 2 ) 2 \k-h\ 2 



Ak 2 



In- 



16fc 3 fci |fe + fcx 



^(Ai), (D5) 



fc 2 (fc 2 -3fc 2 ) fc(fc 2 -fc 2 ) 2 \k-ki\ 



1 1 ^ 



8kf 



1 fc 2 ' 
1 + 



In 



32fcf |fc + fci| 2 



d 3 ki 



rfc 2 + fc 2 (fc 2 -*; 2 ) 2 |fc-fci| 2 
In 



8fc 2 



32fcfc 3 |fc + fci| 2 



12 



Il isl (k)= I d 3 k 1 /3(k 1 ,k-k 1 )/3(-k 1 ,k)P. i ( ° ) (fc 1 ) = --i-fc 2 / d^x-^^i). 



k\ 11 



(D6) 
(D7) 

(D8) 

Vf{k x ), (D9) 
(D10) 



However, only five of the above integrals are independent (for any given i and j), and we have the following relations 

4' S3 (k) = \i 5 i S3 (k) + ±ii' S3 (k), if s \k) = - 1 -i s ; s \k), 4 s \k) = -\4 tS3 {k), ii si {k) = hf si {k). (dii) 

Then we may choose any non-degenerate linear combination of e.g. {/f 5 (fc), . . . , if gJ (fc)} as our set of independent 
integrals. In order to facilitate comparison with the known expressions for the case of a single component fluid [27| . 
we define: 



f[ 5 \k) 



42 



2 

21" 



d 3 kx 

504fc 3 fcf 



3 Ifc — fci I 2 

6fc 7 fci - 79fc 5 fc? + 50fc 3 fc? - 21/cfc 7 + -(fc 2 - fc?) 3 (2fc 2 + 7kf) In ^ 

4 V fc + fci 2 



5 if ai (fc) + y/ 2 ' 51 (*) - (fc) + ^rtf * (fc) + lif* (fc ) 



14 



1 

14" 



1 

14" 



/l* 53 (fc) 



d 3 kx 

168fc 3 fcf 
/f 53 ( 

d 3 k! 



6fc 7 fci - 41fc 5 fc? + 2fc 3 fcf - 3fcfc 7 + ^(fc 2 - fc?) 3 (2fc 2 + fc?) In — 

4 ' |fe + fcx| 



^ } (fci 



( fc) - \i( " (fc) + |if " (fc) - (fc) + 2ir ( fc) 



ifc-fci 



24fc 3 fcf 

d 3 k! 



6fc 7 fci + fc 5 fc 3 + 9fcfc 1 7 + -(fc 2 -fc 2 ) 2 (2fc 4 + 5fc 2 fc 2 +3fc 2 )ln II , , 
4 V |fc + fci| 2 



^(fci) 



^/f 53 (fc) + ^jf &i (k) - u™{k) + °-i s : s3 (k) ^ s3 (k) 



72fc 3 fcf 



6fc 7 fci + 29fc 5 fc? - 18fc 3 fcf + 27fcfc 7 + hk 2 - fc 2 ) 2 (2fc 4 + 9fc 2 fc? + 9k 2 ) In ^ — ^4 

4 |fc + fci| 2 



/f 53 (fc) =ir j (k) 



S l 5 3 i 



(D12) 
(D13) 
(D14) 

(D15) 

^ } (fci) 
(D16) 



-d 3 k!, 



^ 0) (^) : 



as the set of independent integrals appearing in Eq. (|T9"|) . Clearly the functions {/i, /2, /3, /4} are just the {/, g, h, i} 
functions of (27| . while /g is simply a constant. This constant is formally divergent for initial power spectra that 
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fall no faster than fc~ 3 at large fc. However, we regard this to be an artifact of using unsmoothed (i.e. infinitely 
"choppy") initial fields to arbitrarily low spatial scales. Clearly this is unphysical. Hence, in obtaining the numerical 
results presented in this paper, we applied a Gaussian smoothing to all initial power spectra, with smoothing scale 
Smooth = 0.1 h- 1 Mpc, i.e. we set v\f(k) -> v\f {k)e~ k2 . 

Finally, for the sake of completeness, we note that the f(k) and g(k) functions of §VII Al may be expressed with 
the basis integrals as 

f(k) = w 2 ff s °(k) + 2w 1 w 2 ff sb (k)+w 2 2 ff sb (k), (D17) 
g(k) = wlfP\k) + 2w 1 w 2 fP\k) + wlfi' s \k). (D18) 



2. Momentum integrals arising in computing 5P^ MG (k, rf) 

When computing the mode-coupling piece of the one-loop power spectrum, SPj^ MC (k, rj), we encounter the mo- 
mentum integral 

J d 3 k x 7&(k!,k - kx^-ki.k! - k)v£Hh)V$(\k - k x |) . (D19) 

Since a{— ki, — k 2 ) = a(ki, k 2 ) and /?(— kx, — k 2 ) = /?(ki, k 2 ), we see that only six distinct types of integrals can arise 
in the computation. These read 



Jf 5W (fc) = /d 3 k x Wk^k-kOl^^^i^dk-k!!), 



d 3 kl ^!p(0) (fcl) p(0) ( J k ^2kk lX + kl ) , (D2U) 



fc 2 « 

Jf 535fc5 '(fc) = | d 3 k 1 a(k 1 ,k-k 1 )a(k-k 1 ,k 1 )^ 0) (fc 1 )^° ) (|k-k 1 |) , 

k 2 x(k-kxx) r>(a) ru ^(0) 
jfex (fc 2 - 2/fc/cia; + fc 2 ) 



= / d J k x - ^ L— P%>( kl )r%> ( \jk 2 - 2fcfc^ + fc 2 ) , (D21) 



= / d 'ki „ „„ v l — ' , 2 , 2 n; ; (fei)nr ( \/fc 2 - 2*fci* + fc 2 ) , (D23) 



jf 4 ^ 4 '^) = y'd 3 k 1 a(k 1 ,k-k 1 ) / 3(k 1 ,k-k 1 )7'§ ) (fei)^ ) (|k-k 1 |) , 

- / ^ 2ki^i^ x \k^ ^ (v*-^+*0 ■ ^ 

k 4 (h-kx) 2 „ i0)fu ^(o) 
4fc 2 {k 2 ~ 2kk lX + k'f) 

Jr 5kS '(k) = |d 3 k! Kk-k^kOl'^^i^dk-kil), 

k 2 (k-kix) 2 ^(Q)^ ^(0) 

(fc 2 - 2fcfci2i + A: 2 ) 

jf«^ fc * ! (fc) = /"^kiaCk-ki.ki^fki.k-ki^ffei^dk-k! 

fc 3 [zfc 2 - fci (x 2 + 1) fc + fc 2 x] ^,(0)^ ^(0) 
2fci (fc 2 - 2fcfcix + fc 2 ) 



d 3 kl - b-f-^ ^fe)^ V fc2 - 2fc ^ + fc2 ) ' < D24 ) 



However, not all of these integrals are independent, since changing the variable of integration as ki — > k' x = k — ki 
merely exchanges jf 43W (fc) with jf 4W (fc) and jf 4 ' 4 * 4 ' (fc) with jf 4W (fc). Thus, for (tj) = («), ^5 and J 6 
are not independent integrals: 



J 4W43 (fc) = jf 4W (fc), J 4!434 ' 43 (fc) = jf 434,43 (fc). (D26) 



3G 



jf ^ 3 ( k ) = jf * w (k) , jf s ' 6i5i (k) = 4 S3&k&l (k) , 4 k&l&lS3 (k) = 4 &3Sk&l (k) 
jfs'^ik) = 4' S36ks '(k) , 4 Wsi (k) = j[ S3Skg '(k) . 4 k5 ' 5 * S3 (k) = 4 i5i5ks '{k) . 



rS i 8 j 8 i 8 j 
In 


(k) 




n = l,.. 


..,4, 


(ii) = {(cc),(cb),(bb)}, 


r8 i 8 j 8 k 8 l 
In 


(ft) 


= 4 P5k5l (k), 


n = 1, . 


..,6, 


= {(cc)(cb), (cc)(bb), (cb)(bb)} 



On the other hand, if (ij) 7^ (fcZ), then we find that we may express each J* 5 s &] (fc) with some s (fc): 

(D27) 

°- (fc) = J° °" ° " (fc) , J 5 ° ° ° (fc) = J£ °" " " (fc) . Jg° ° ° °" (fc) = J 3 ° " " (fc) . 

Thus, we may define 

(D28) 
(D29) 

as the set of independent integrals appearing in Eq. ([85f . 

Finally, for the sake of completeness, we note that the F(k), G(k) and if (fc) functions of WIIBI mav be expressed 
with these basis integrals as 

F{k) = W6 J ' m{k) + 1 ffm(k) + S + J| ^ (fc) + W6 J ' m{k) + ^ m{k) ' (D30) 

G(fc) = lie Jl ^ (fc) + ^ J * m W + ^4 m (k) + l pl m {k) + j^4 m (k) + H Jp(fc) , (D31) 
* (*) = + + § + ^4 m {k) + f 8 4 m (k) + I Jp(fc) , (D32) 

where 

4««(jfe) = w fjf 5W (fc) + 2™<Wf ^"(fc) + ^Jf wl ' b (t) 

+ 2^ 2 Jf ^(fc) + 4 W > 2 2 jf ^'(fc) + 2 Wl ^ 2 3 Jf 5W (fc) (D33) 
+ w\wl J r f sbgcsc (k) + 2^1 jf 5W (fc) + wfjf * W (fc) . 

In Eqs (|D30| - (jD32)) above, we clearly have 4 m (k) = jf sss (fc) and J<p s (fc) = jf sss (fc). Ho wever , to ob tain the 
compact expressions for F(fc), G(fc) and -ff(fc) as in Eqs (fi?5|) - (|100p , one must substitute Eqs (|D20p - (|D25|) a s they 
stand into the general formulae, Eqs (|D30|) - (|D32[) . disregarding the above identities. 

APPENDIX E: THE LINEAR PROPAGATOR FOR N = 3 

For N = 3, the propagator takes the form as in Eq. (f2"T|) : 

9ab(v) ^J2 elV 9ab,l, (El) 
I 
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with I e {1,0, —1/2, —3/2}. The g a b,i are constant matrices: 
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9ab,-3/2 



The eigenvectors of the linear propagator g a b are: 
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-2wi 
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-2wi 
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3wi 
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-3w 3 


3w 3 _ 



,(1) 



1 
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1 
1 

V i/ 



,(2) 
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-1 
2/3 
-1 
2/3 

V -i / 



,(3,1) 



/ W 2 \ 


-Wl 




V o / 



,(3,2) 



/ w 3 \ 




-Wl 

V o J 
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-2wi 
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The eigenvectors correspond to the following eigenvalues: 



u 



« : Ai 



u 



< 2 ) : A 2 = e' 3 "/ 2 , 



u 



< 3 < m ) : A 3 = 1 , : A 4 = ■ 



-n/2 



(E2) 



(E3) 



(E4) 



(E5) 



,(4,2) _ 



/ 2w 3 \ 

-w 3 



-2wi 
\ Wl / 
(E6) 



m = 1,2. 



(E7) 



APPENDIX F: THE LINEAR PROPAGATOR FOR GENERAL N 

Here we compute the linear propagator, g a b(i]) = [&ab(s), s, 77], for general N, where £ _1 [/(s), s, 77] denotes the 
inverse Laplace transform of the function f(s) from the variable s to the variable 77 and a a b(s) = [s5 a b + Qab)^ 1 - 
We begin by noting that 



(si + n)(2j-l)(2fc-l) = sdjk , 
(si + n)(2j-l)(2fe) = -$jk j 



- 2 w k , 



(si + n)(2j)(2fc-i) 

(si + ^)(2j)(2fe) = Q + , 



(Fl) 

(F2) 
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with j,k — 1,2, ... ,N. Now we claim that a a b{s) may be written in the following form for general N: 



3w k 



S 



jk 



* W -l)(2fc-l)(*) - ^_ 1)(2s + 3) ' s > 

3wk 



0"(2j)(2fe-l)(s) 
cr (2j-l)(2fe)(s) 
C r (2j)(2fc)(s) 



(s-l)(2s + 3) ' 

6wfe 

s(s-l)(2s + l)(2s + 3) ' s(2s + l) 



2(5 



'j'fc 



+ ■ 



2(5 



(s-l)(2s + l)(2s + 3) (2s + 1) ' 
(j, fc = 1, 2, . . . , iV). We prove this claim by explicit computation. We calculate 

2N N 

[(j(s) • (sl + 0)] ab = ^er ac (s)(sl + ft) cfc = [o- (2fe-l)(s)(sl + n)(2fe_l)6 + 0- (2fe)(s)(sl + n)( 2 fe)6] 
c=l fc=l 



(F3) 



(F4) 



where in the last expression, we have split the sum over c into two sums, corresponding to odd and even c respectively. 
Now we examine the four cases corresponding to a and b being odd or even separately. We find 

1. a = In- 1, b = 2m- 1 

AT 

^ [ cr (2n-l)(2fc-l)( s )( s l + ^)(2fe-l)(2m-l) + cr (2n-l)(2fe) ( s ) (si + ^)(2fe)(2m-l)] 



fc=l 
N 

=E 

fc=i 



3w fc (5„fc , 

H (sdfcm) + 



s(s-l)(2s + 3) s 



+ 



3w m 3 

(s - l)(2s + 3) 2 



s(s-l)(2s + l)(2s + 3) s(2s + l) 



+ 



-2^ 



s(s-l)(2s + l)(2s + 3) s(2s + l) 



(F5) 



where we have used that X^fe=i w fe = 1- 
2. a = 2n, 6 = 2m- 1 



^ [c(2n)(2fc-l)(s)(sl + ^)(2fe-l)(2m-l) + f(2n)(2fe) (s)(sl + ^)(2fc)(2m-l)] 

2(5„fc 



fc=i 

N 

E 



3w fc 
(«-l)(2* + 3) 

3sw m 3 



(s&m) 



(s - l)(2s + l)(2s + 3) + (2s + 1) 
6 2 



(s-l)(2s + 3) 2 V(s - l)(2s + l)(2s + 3) (2s + 1) 
=0, 



where we have used that Y^k=i Wk = ^■ 
3. a = 2n - 1, b = 2m 



3 

-2 W " 



(F6) 



AT 



^ [f(2n-l)(2fe-l)(s)(sl + 0)(2fe-l)(2m) + °"(2n-l)(2fc) (s){sl + fi)(2fc)(2ro)] 



fe=l 
AT 



3wfc + ^ ] (-4™) + 



s(s-l)(2s + 3) s 
3w m 5. 



6wi 



26 



n k 



s(s-l)(2s + 3) s 



s(s-l)(2s + l)(2s + 3) s(2s + l) 
Qw m 2,d nm 



s(s-l)(2s + l)(2s + 3) s(2s + l) 



=0. 



- + s ) 4™ 



(F7) 
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4. a = 2n, b = 2m 

N 

^ [f(2n)(2fc-l)(s)(sl + fi)(2fc-l)(2m) + C(2n)(2fc) + fi)(2fc)(2m)] 

fc=l 
iV 

=E 

fe=i 

3w m /l \ / 6w m 2S nm \ 

(s-l)(2a + 3) \2 7 V(s-l)(2s + l)(2s + 3) (2s + 1) J 

Finally, we recognize that performing the inverse Laplace transform of Eq. (|F3|). we simply obtain Eq. 
Similar explicit computations (which however we do not reproduce) establish that for general N, the propagator 
has the following 27V eigenvalues and eigenvectors: 

1. Eigenvalue Ai = e 1 ' with multiplicity one. The corresponding eigenvector is: 

nW = (l,l,...,l,l) T . (F9) 

2. Eigenvalue A2 = e~ 3r) ^ 2 with multiplicity one. The corresponding eigenvector is: 

= (2/3,-1,..., 2/3, -1) T . (F10) 

3. Eigenvalue A3 = 1 with multiplicity N — 1. The corresponding eigenvectors are: 

u (3,m) _ ( Wmj o, ... ,0, -wi,0, . . . ,0) T = w m S al - wi5 a(2m -i) , (Fll) 

(rn = 2, ... , AT), i.e. the two non-zero entries in u^ 3,m - ) are the first and the 2to — 1-th. 

4. Eigenvalue A3 = e~ v ^ 2 with multiplicity N — 1. The corresponding eigenvectors are: 

u (4, m ) = (2 Wmj -iy m , 0, . . . , 0, -2wi, wx, 0, . . . , 0) T = 2w m (5 Q i - w m 5 a2 - 2wx5 a ^ m -x) + wiS a{2m) , (F12) 
(to = 2, ... , AT), i.e. the four non-zero entries in u^ 3,m ^ are the first and second and the 2to — 1-th and 2m-th. 
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(F8) 



